gestural_alignment_rate
Last updated: 2025-11-24
library(tidyverse)
library(plotrix) #std.error
library(hypr) #for generating contrast coding
library(ggthemes)
library(ggh4x)
library(dagitty) #dagitty for causal diagrams
library(DiagrammeR) #grViz for causal diagrams
library(DiagrammeRsvg)
library(rsvg)
library(plotly) #ggplotly for interactive plots
library(brms) #bayesian regression
library(bayestestR) #p_direction, hdi
library(tidybayes) #add_epred_draws
library(emmeans)
library(doParallel) #set the number of cores
library(ppcor) #partial correlation
library(svglite)
### Set global options
options(digits = 3) # set the default number of digits to 3
cores = as.integer(detectCores(logical = FALSE) * 0.8) # set the number of cores to use
registerDoParallel(cores=cores) # register the number of cores to use for parallel processing
options(mc.cores = cores)
options(brms.backend = "cmdstanr") #this will speed up the model fitting
### MCMC options
niter = 20000 #number of iterations
nwu = 2000 #number of warmups
### Rmd settings
knitr::opts_chunk$set(echo=TRUE, warning=FALSE, message=FALSE, fig.path="figures_md/speakerB/gest_alignment/")### Custom functions
########### pp_check_each_round ############
pp_check_each_round <- function(m, data, round_info) {
df_sub = filter(data, round == round_info)
p = pp_check(m,
type = "bars",
ndraws = 100,
newdata = df_sub) +
coord_cartesian(xlim = c(0, 10)) +
ggtitle(round_info)
return(p)
}
########### pp_check_each_condition ############
pp_check_each_condition <- function(m, data, condition_info) {
df_sub = filter(data, condition == condition_info)
p = pp_check(m,
type = "bars",
ndraws = 100,
newdata = df_sub) +
coord_cartesian(xlim = c(0, 10)) +
ggtitle(condition_info)
return(p)
}
########### plot_posterior ############
plot_posterior <- function(model, model2=NA, model3=NA,
interaction=FALSE, include_intercept=FALSE,
xlim_cond=1.5, xlim_round=2, xlim_lex=0.4){
### extract the posterior draws
posterior_beta1 <- model %>%
gather_draws(`b_.*`, regex = TRUE) %>%
mutate(intercept = str_detect(.variable, "Intercept"),
component = ifelse(str_detect(.variable, ":"), "Interaction",
ifelse(str_detect(.variable, "round"), "Round",
ifelse(str_detect(.variable, "Intercept"), "Intercept",
ifelse(str_detect(.variable, "lex_align"), "N. lex align",
"Visibility")))))
if (length(model2) == 1){ #if model2 is NA
posterior_beta = posterior_beta1
} else {
posterior_beta1 <- posterior_beta1 %>%
filter(.variable != "b_lex_align_c")
posterior_beta2 <- model2 %>%
gather_draws(`b_.*`, regex = TRUE) %>%
filter(.variable == "b_lex_align_c") %>%
mutate(component = "N. lex align")
posterior_beta3 <- model3 %>%
gather_draws(`b_.*`, regex = TRUE) %>%
filter(.variable == "b_condition_sumAO_Sym") %>%
mutate(component = "Visibility")
posterior_beta = rbind(posterior_beta1, posterior_beta2, posterior_beta3)
}
posterior_beta = posterior_beta %>%
mutate(.variable = recode(.variable,
"b_Intercept" = "Intercept",
"b_conditionAsymAV" = "SymAV--AsymAV",
"b_conditionAO" = "SymAV--AO",
"b_conditionAsym_Sym" = "AsymAV--SymAV",
"b_conditionAO_Asym" = "AO--AsymAV",
"b_condition_sumAO_Sym" = "AO--SymAV",
"b_lex_align_c" = "N. lex align",
"b_round_c" = "Centered round",
"b_log_round_c" = "Centered log(round)",
"b_roundR12" = "R1--R2",
"b_roundR23" = "R2--R3",
"b_roundR34" = "R3--R4",
"b_roundR45" = "R4--R5",
"b_roundR56" = "R5--R6",
"b_conditionAsymAV:round1" = "SymAV--AsymAV: R1--R2",
"b_conditionAsymAV:round2" = "SymAV--AsymAV: R2--R3",
"b_conditionAsymAV:round3" = "SymAV--AsymAV: R3--R4",
"b_conditionAsymAV:round4" = "SymAV--AsymAV: R4--R5",
"b_conditionAsymAV:round5" = "SymAV--AsymAV: R5--R6",
"b_conditionAO:round1" = "SymAV--AO: R1--R2",
"b_conditionAO:round2" = "SymAV--AO: R2--R3",
"b_conditionAO:round3" = "SymAV--AO: R3--R4",
"b_conditionAO:round4" = "SymAV--AO: R4--R5",
"b_conditionAO:round5" = "SymAV--AO: R5--R6",
"b_conditionAsym_Sym:round_c" = "Centered round: Asym--Sym",
"b_conditionAO_Sym:round_c" = "Centered round: AO--Sym",
"b_conditionAO_Asym:round_c" = "Centered round: AO--Asym",
"b_conditionAsym_Sym:log_round_c" = "Centered log(round): Asym--Sym",
"b_conditionAO_Sym:log_round_c" = "Centered log(round): AO--Sym",
"b_conditionAO_Asym:log_round_c" = "Centered log(round): AO--Asym"),
.variable = factor(.variable,
levels = c("AO--SymAV", "AO--AsymAV", "AsymAV--SymAV",
"R1--R2", "R2--R3", "R3--R4", "R4--R5", "R5--R6",
"N. lex align")),
component = factor(component,
levels = c("Intercept", "Visibility", "Round",
"Interaction", "N. lex align")))
if (include_intercept == F){
posterior_beta = posterior_beta %>% filter(component != "Intercept")
}
### change variables if only main effects are plotted
if (interaction == F) {
posterior_beta = filter(posterior_beta, !str_detect(.variable, ":"))
fill_manual_values = c("steelblue", "steelblue", "steelblue")
} else{
fill_manual_values = c("steelblue", "steelblue", "steelblue", "steelblue")
}
### plot the posterior distributions
p_posterior = ggplot(posterior_beta,
aes(x = .value, y = fct_rev(.variable),
fill = component)) +
geom_vline(xintercept = 0, size = 1) +
stat_halfeye(aes(slab_alpha = intercept),
normalize = "panels",
slab_alpha = .5,
.width = c(0.89, 0.95),
point_interval = "median_hdi") +
scale_fill_manual(values = fill_manual_values) +
scale_slab_alpha_discrete(range = c(0.8, 0.4)) +
guides(fill = "none", slab_alpha = "none") +
labs(x = "Coefficient", y = "Effect") +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "none",
strip.text = element_text(size = 15, face = 'bold'),
strip.background = element_blank(),
panel.grid.major.x = element_line(color = "grey90",
linetype = "solid",
size = 0.5),
panel.grid.major.y = element_line(color = "grey90",
linetype = "solid",
size = 0.5),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines")) +
facet_wrap(vars(component), ncol = 3, scales = "free") +
facetted_pos_scales(
x = list(
scale_x_continuous(limits = c(-xlim_cond, xlim_cond)),
scale_x_continuous(limits = c(-xlim_round, xlim_round)),
scale_x_continuous(limits = c(-xlim_lex, xlim_lex))))
return(p_posterior)
}
########### pp_update_plot ############
pp_update_plot <- function(post_sample, model_type="nb", interaction=TRUE){
sum = ifelse("b_condition_sumAO_Sym" %in% colnames(post_sample), T, F)
round_bd = ifelse("b_roundR12" %in% colnames(post_sample), T, F)
intercept = ggplot(post_sample) +
geom_density(aes(prior_Intercept), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_Intercept), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Intercept') +
theme_classic()
### Visibility condition
if (sum == F){
cond1 = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_conditionAsym_Sym), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Asym--Sym') +
theme_classic()
cond2 = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_conditionAO_Asym), fill="#FC4E07", color="black",alpha=0.6) +
xlab('AO--Asym') +
theme_classic()
} else {
cond1 = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_condition_sumAO_Sym), fill="#FC4E07", color="black",alpha=0.6) +
xlab('AO--Sym') +
theme_classic()
cond2 = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_condition_sumAsym_Sym), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Asym--Sym') +
theme_classic()
}
### Round
if (interaction) {
if (round_bd){
r1 = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_roundR12), fill="#FC4E07", color="black",alpha=0.6) +
xlab('R1--R2') +
theme_classic()
r2 = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_roundR23), fill="#FC4E07", color="black",alpha=0.6) +
xlab('R2--R3') +
theme_classic()
r3 = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_roundR34), fill="#FC4E07", color="black",alpha=0.6) +
xlab('R3--R4') +
theme_classic()
r4 = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_roundR45), fill="#FC4E07", color="black",alpha=0.6) +
xlab('R4--R5') +
theme_classic()
r5 = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_roundR56), fill="#FC4E07", color="black",alpha=0.6) +
xlab('R5--R6') +
theme_classic()
} else {
if ("b_round_c" %in% colnames(post_sample)) {
round = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_round_c), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Centered round') +
theme_classic()
if (sum == F){
cond1_round = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(`b_conditionAsym_Sym:round_c`), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Centered Round: Asym--Sym') +
theme_classic()
cond2_round = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(`b_conditionAO_Asym:round_c`), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Centered Round: AO--Asym') +
theme_classic()
} else {
cond1_round = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(`b_condition_sumAO_Sym:round_c`), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Centered Round: AO--Sym') +
theme_classic()
cond2_round = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(`b_condition_sumAsym_Sym:round_c`), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Centered Round: Asym--Sym') +
theme_classic()
}
} else {
round = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_log_round_c), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Centered log(round)') +
theme_classic()
if (sum == F){
cond1_round = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(`b_conditionAsym_Sym:log_round_c`), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Centered log(round): Asym--Sym') +
theme_classic()
cond2_round = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(`b_conditionAO_Asym:log_round_c`), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Centered log(round): AO--Asym') +
theme_classic()
} else {
cond1_round = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(`b_condition_sumAO_Sym:log_round_c`), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Centered log(round): AO--Sym') +
theme_classic()
cond2_round = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(`b_condition_sumAsym_Sym:log_round_c`), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Centered log(round): Asym--Sym') +
theme_classic()
}
}
}
}
if (model_type == "nb"){
shape = ggplot(post_sample) +
geom_density(aes(prior_shape), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(shape), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Shape') +
scale_x_continuous(limits = c(0, 10)) +
theme_classic()}
else if (model_type == "zinb") {
shape = ggplot(post_sample) +
geom_density(aes(prior_shape), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(shape), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Shape') +
scale_x_continuous(limits = c(0, 10)) +
theme_classic()
zi = ggplot(post_sample) +
geom_density(aes(prior_zi), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(zi), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Zero-inflation') +
theme_classic()}
else if (model_type == "zibt"){
phi = ggplot(post_sample) +
geom_density(aes(prior_phi), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(phi), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Precision') +
theme_classic()
zoi = ggplot(post_sample) +
geom_density(aes(prior_zoi), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(zoi), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Zero-inflation') +
theme_classic()
coi = ggplot(post_sample) +
geom_density(aes(prior_coi), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(coi), fill="#FC4E07", color="black",alpha=0.6) +
xlab('One-inflation') +
theme_classic()
}
### display the plots
if (interaction==F){
if (model_type=="nb"){
gridExtra::grid.arrange(intercept, cond1, cond2, shape, ncol=2)}
else if (model_type=="zinb"){
gridExtra::grid.arrange(intercept, cond1, cond2, shape, zi, ncol=3)}
else if (model_type=="zibt"){
gridExtra::grid.arrange(intercept, cond1, cond2, phi, zoi, coi, ncol=3)
}
} else {
if (round_bd == T){
if (model_type=="nb"){
gridExtra::grid.arrange(intercept, cond1, cond2, shape, ncol=3)}
else if (model_type=="zinb"){
gridExtra::grid.arrange(intercept, cond1, cond2, shape, zi, ncol=3)}
else if (model_type=="zibt"){
gridExtra::grid.arrange(intercept, cond1, cond2, phi, zoi, coi, ncol=3)
}
} else {
if (model_type=="nb"){
gridExtra::grid.arrange(intercept, cond1, cond2, round,
cond1_round, cond2_round, shape, ncol=3)}
else if (model_type=="zinb"){
gridExtra::grid.arrange(intercept, cond1, cond2, round,
cond1_round, cond2_round, shape, zi, ncol=3)}
else if (model_type=="zibt"){
gridExtra::grid.arrange(intercept, cond1, cond2, round,
cond1_round, cond2_round, phi, zoi, coi, ncol=3)
}
}
}
}1. Data preparation
1.1 Load data
### trial_info_speaker.csv
df_trial_info = read.csv("data/trial_info_speaker.csv", stringsAsFactors = TRUE) %>%
filter(num_words != 0) %>% # remove trials that were accidentally skipped)
mutate(pair = factor(pair),
role = factor(ifelse(speaker == director, "director", "matcher")),
round_c = as.integer(round) - mean(as.integer(round)),
round_n = factor(round),
round = factor(paste0("R", round)),
log_round = log(as.integer(round)),
log_round_c = log_round - mean(log_round),
condition = factor(condition,
levels = c("Sym", "Asym", "AO"),
labels = c("SymAV", "AsymAV", "AO")),
condition_sum = condition,
duration_s = duration_ms / 1000,
n_words_100 = num_words / 100,
n_iconic_c = num_iconic_gestures - mean(num_iconic_gestures),
n_iconic_per_100words = num_iconic_gestures / n_words_100,
n_iconic_100 = num_iconic_gestures / 100,
lex_align_c = num_lex_align - mean(num_lex_align),
lex_align_rate = num_lex_align / n_words_100,
gest_align_rate_100words = num_gestural_align / n_words_100,
gest_align_rate = num_gestural_align / num_iconic_gestures) %>%
rename(trial_duration_s = duration_s, trial_duration_m = duration_m) %>%
filter(speaker == "B")
### condition info
df_condition_info = read.csv("data/condition_info.csv", stringsAsFactors = TRUE) %>%
mutate(pair = factor(pair),
condition = factor(condition,
levels = c("Sym", "Asym", "AO"),
labels = c("SymAV", "AsymAV", "AO")))1.2 Summarize data
summarize_demographics <- function(df) {
df %>%
summarize(total_s = sum(trial_duration_s),
total_m = total_s / 60,
### trial duration ###
mean_trial_dur_s = mean(trial_duration_s),
mean_trial_dur_m = mean(trial_duration_m),
sd_trial_dur_m = sd(trial_duration_m),
lci_trial_dur_m = mean_trial_dur_m - 1.96 * sd_trial_dur_m / sqrt(n()),
uci_trial_dur_m = mean_trial_dur_m + 1.96 * sd_trial_dur_m / sqrt(n()),
### words ###
# number of words
num_words_total = sum(num_words),
mean_num_words = mean(num_words),
num_words_100 = mean(num_words / 100),
sd_num_words = sd(num_words),
lci_num_words = mean_num_words - 1.96 * sd_num_words / sqrt(n()),
uci_num_words = mean_num_words + 1.96 * sd_num_words / sqrt(n()),
num_words_per_min = num_words_total / total_m,
# number of content words
num_content_total = sum(num_content_words),
mean_num_content = mean(num_content_words),
sd_num_content = sd(num_content_words),
lci_num_content = mean_num_content - 1.96 * sd_num_content / sqrt(n()),
uci_num_content = mean_num_content + 1.96 * sd_num_content / sqrt(n()),
num_content_per_min = num_content_total / total_m,
### lexical alignment ###
# raw frequency
num_lex_align_total = sum(num_lex_align, na.rm = T),
mean_num_lex_align = mean(num_lex_align, na.rm = T),
sd_num_lex_align = sd(num_lex_align, na.rm = T),
lci_num_lex_align = mean_num_lex_align - 1.96 * sd_num_lex_align / sqrt(n()),
uci_num_lex_align = mean_num_lex_align + 1.96 * sd_num_lex_align / sqrt(n()),
num_lex_align_per_min = num_lex_align_total / total_m,
# rate per 100 words
mean_lex_align_per_100words = mean(lex_align_rate, na.rm=T),
sd_lex_align_per_100words = sd(lex_align_rate, na.rm=T),
se_lex_align_per_100words = std.error(lex_align_rate, na.rm=T),
lci_lex_align_per_100words = mean_lex_align_per_100words - 1.96 * se_lex_align_per_100words,
uci_lex_align_per_100words = mean_lex_align_per_100words + 1.96 * se_lex_align_per_100words,
### iconic gestures ###
# raw frequency
num_iconic_total = sum(num_iconic_gestures, na.rm = T),
mean_num_iconic = mean(num_iconic_gestures, na.rm = T),
sd_num_iconic = sd(num_iconic_gestures, na.rm = T),
lci_num_iconic = mean_num_iconic - 1.96 * sd_num_iconic / sqrt(n()),
uci_num_iconic = mean_num_iconic + 1.96 * sd_num_iconic / sqrt(n()),
num_iconic_per_min = num_iconic_total / total_m,
# rate per 100 words
mean_iconic_per_100words = mean(n_iconic_per_100words, na.rm=T),
sd_iconic_per_100words = sd(n_iconic_per_100words, na.rm=T),
se_iconic_per_100words = std.error(n_iconic_per_100words, na.rm=T),
lci_iconic_per_100words = mean_iconic_per_100words - 1.96 * se_iconic_per_100words,
uci_iconic_per_100words = mean_iconic_per_100words + 1.96 * se_iconic_per_100words,
### gestural alignment ###
# raw frequency
num_gest_align_total = sum(num_gestural_align, na.rm = T),
mean_num_gest_align = mean(num_gestural_align, na.rm = T),
sd_num_gest_align = sd(num_gestural_align, na.rm = T),
lci_num_gest_align = mean_num_gest_align - 1.96 * sd_num_gest_align / sqrt(n()),
uci_num_gest_align = mean_num_gest_align + 1.96 * sd_num_gest_align / sqrt(n()),
num_gest_align_per_min = num_gest_align_total / total_m,
# rate per 100 words
mean_gest_align_per_100words = mean(gest_align_rate_100words, na.rm=T),
sd_gest_align_per_100words = sd(gest_align_rate_100words, na.rm=T),
se_gest_align_per_100words = std.error(gest_align_rate_100words, na.rm=T),
lci_gest_align_per_100words = mean_gest_align_per_100words - 1.96 * se_gest_align_per_100words,
uci_gest_align_per_100words = mean_gest_align_per_100words + 1.96 * se_gest_align_per_100words,
# rate per iconic gestures
mean_gest_align_prop = mean(gest_align_rate, na.rm=T),
sd_gest_align_prop = sd(gest_align_rate, na.rm=T),
se_gest_align_prop = std.error(gest_align_rate, na.rm=T),
lci_gest_align_prop = mean_gest_align_prop - 1.96 * se_gest_align_prop,
uci_gest_align_prop = mean_gest_align_prop + 1.96 * se_gest_align_prop,
### number of trials ###
trial_n = n()) %>%
ungroup()
}
summarize_dur <- function(df){
df %>%
summarize(mean_dur_m = mean(total_m),
sd_dur_m = sd(total_m),
se_dur_m = std.error(total_m),
lci_dur_m = mean_dur_m - 1.96 * se_dur_m,
uci_dur_m = mean_dur_m + 1.96 * se_dur_m) %>%
ungroup()
}mean by condition
### demographics by pair
trial_info_pair = df_trial_info %>%
group_by(pair) %>%
summarize_demographics() %>%
left_join(., df_condition_info) %>%
dplyr::select(pair, condition, total_s:trial_n)
### calculate mean experiment duration
mean_dur_cond = trial_info_pair %>%
group_by(condition) %>%
summarize_dur()
### summary statistics
trial_info_cond = df_trial_info %>%
group_by(condition) %>%
summarize_demographics() %>%
left_join(., mean_dur_cond) %>%
dplyr::select(condition, total_m, mean_dur_m:uci_dur_m,
everything(),
-ends_with("_s"))
trial_info_cond### export it as csv
tbl = trial_info_cond %>%
dplyr::select(!starts_with(c("sd_", "se_", "num_", "trial_"))
& !contains("total")) %>%
pivot_longer(cols = !condition,
names_to = "name",
values_to = "value") %>%
mutate(stats = sub("_.*", "", name),
name = sub("[[:alpha:]]+_", "", name),
across(where(is.numeric), ~ str_squish(format(round(., 2), scientific = FALSE)))) %>%
pivot_wider(names_from = c("condition", "stats"),
values_from = "value") %>%
mutate(SymAV_ci = paste0("[", SymAV_lci, ", ", SymAV_uci, "]"),
AsymAV_ci = paste0("[", AsymAV_lci, ", ", AsymAV_uci, "]"),
AO_ci = paste0("[", AO_lci, ", ", AO_uci, "]")) %>%
dplyr::select(name, SymAV_mean, SymAV_ci, AsymAV_mean, AsymAV_ci, AO_mean, AO_ci)
write_csv(tbl, "tables/descriptive_B.csv")mean by round
trial_info_pair_round = df_trial_info %>%
group_by(pair, round, round_n) %>%
summarize_demographics() %>%
left_join(., df_condition_info)
### calculate mean round duration
mean_dur_round = trial_info_pair_round %>%
group_by(round, round_n) %>%
summarize_dur()
trial_info_round = df_trial_info %>%
group_by(round, round_n) %>%
summarize_demographics() %>%
left_join(., mean_dur_round) %>%
dplyr::select(round, round_n, total_m, mean_dur_m:uci_dur_m,
everything(),
-ends_with("_s"))
trial_info_roundmean by condition x round
### calculate mean duration by condition x round
mean_dur_cond_round = trial_info_pair_round %>%
group_by(condition, round, round_n) %>%
summarize_dur()
trial_info_cond_round = df_trial_info %>%
group_by(condition, round, round_n) %>%
summarize_demographics() %>%
left_join(., mean_dur_cond_round) %>%
dplyr::select(condition, round, round_n,
total_m, mean_dur_m:uci_dur_m,
everything(),
-ends_with("_s"))
trial_info_cond_round3. Contrast coding
### visibility condition: difference coding
h_cond = hypr(AO_Asym = AsymAV ~ AO,
Asym_Sym = SymAV ~ AsymAV,
levels = levels(df_trial_info$condition))
h_cond## hypr object containing 2 null hypotheses:
## H0.AO_Asym: 0 = AsymAV - AO
## H0.Asym_Sym: 0 = SymAV - AsymAV
##
## Call:
## hypr(AO_Asym = ~AsymAV - AO, Asym_Sym = ~SymAV - AsymAV, levels = c("SymAV",
## "AsymAV", "AO"))
##
## Hypothesis matrix (transposed):
## AO_Asym Asym_Sym
## SymAV 0 1
## AsymAV 1 -1
## AO -1 0
##
## Contrast matrix:
## AO_Asym Asym_Sym
## SymAV 1/3 2/3
## AsymAV 1/3 -1/3
## AO -2/3 -1/3contrasts(df_trial_info$condition) = contr.hypothesis(h_cond)
### visibility condition: sum coding
h_cond = hypr(AO_Sym = SymAV ~ AO,
Asym_Sym = SymAV ~ AsymAV,
levels = levels(df_trial_info$condition))
h_cond## hypr object containing 2 null hypotheses:
## H0.AO_Sym: 0 = SymAV - AO
## H0.Asym_Sym: 0 = SymAV - AsymAV
##
## Call:
## hypr(AO_Sym = ~SymAV - AO, Asym_Sym = ~SymAV - AsymAV, levels = c("SymAV",
## "AsymAV", "AO"))
##
## Hypothesis matrix (transposed):
## AO_Sym Asym_Sym
## SymAV 1 1
## AsymAV 0 -1
## AO -1 0
##
## Contrast matrix:
## AO_Sym Asym_Sym
## SymAV 1/3 1/3
## AsymAV 1/3 -2/3
## AO -2/3 1/3contrasts(df_trial_info$condition_sum) = contr.hypothesis(h_cond)
### round
bacward_diff = hypr(R12 = R2 ~ R1,
R23 = R3 ~ R2,
R34 = R4 ~ R3,
R45 = R5 ~ R4,
R56 = R6 ~ R5,
levels = levels(df_trial_info$round))
bacward_diff## hypr object containing 5 null hypotheses:
## H0.R12: 0 = R2 - R1
## H0.R23: 0 = R3 - R2
## H0.R34: 0 = R4 - R3
## H0.R45: 0 = R5 - R4
## H0.R56: 0 = R6 - R5
##
## Call:
## hypr(R12 = ~R2 - R1, R23 = ~R3 - R2, R34 = ~R4 - R3, R45 = ~R5 -
## R4, R56 = ~R6 - R5, levels = c("R1", "R2", "R3", "R4", "R5",
## "R6"))
##
## Hypothesis matrix (transposed):
## R12 R23 R34 R45 R56
## R1 -1 0 0 0 0
## R2 1 -1 0 0 0
## R3 0 1 -1 0 0
## R4 0 0 1 -1 0
## R5 0 0 0 1 -1
## R6 0 0 0 0 1
##
## Contrast matrix:
## R12 R23 R34 R45 R56
## R1 -5/6 -2/3 -1/2 -1/3 -1/6
## R2 1/6 -2/3 -1/2 -1/3 -1/6
## R3 1/6 1/3 -1/2 -1/3 -1/6
## R4 1/6 1/3 1/2 -1/3 -1/6
## R5 1/6 1/3 1/2 2/3 -1/6
## R6 1/6 1/3 1/2 2/3 5/6contrasts(df_trial_info$round) = contr.hypothesis(bacward_diff)
### role (director / matcher)
contrasts(df_trial_info$role) = c(-0.5, 0.5)
contrasts(df_trial_info$role)## [,1]
## director -0.5
## matcher 0.5==== Iconic gestures ====
4. Number of iconic gestures
4.1 DataViz: frequency
bp: total by condition
bp_iconic_by_cond = ggplot(data=trial_info_pair,
aes(x=condition, y=num_iconic_total, fill=condition)) +
geom_jitter(aes(color = pair),
size = 1, alpha = 1,
width = 0.07, height = 0) +
geom_boxplot(width = .2,
outlier.shape = NA, alpha = 0.7) +
scale_fill_manual(values = c("#ED6B06", "#00786A", "darkgrey")) +
labs(x = "Visibility",
y = "Total N of iconic gestures per pair") +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "none",
strip.text = element_text(size = 15, face = 'bold'),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines"))
ggplotly(bp_iconic_by_cond)bp: mean by condition
bp_mean_iconic_by_cond = ggplot(data=trial_info_pair,
aes(x=condition, y=mean_num_iconic, fill=condition)) +
geom_jitter(aes(color = pair),
size = 1, alpha = 1,
width = 0.07, height = 0) +
geom_boxplot(width = .2,
outlier.shape = NA, alpha = 0.7) +
geom_point(data = trial_info_cond,
aes(y = mean_num_iconic),
shape = 21, size = 3, fill = "white") +
scale_y_continuous(limits = c(0, 4), breaks = seq(0, 4, 1)) +
scale_fill_manual(values = c("#ED6B06", "#00786A", "darkgrey")) +
labs(x = "Visibility",
y = "Mean N. iconic gest per trial") +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "none",
strip.text = element_text(size = 15, face = 'bold'),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines"))
ggplotly(bp_mean_iconic_by_cond)bp: mean by condition x round
bp_mean_iconic_by_round_cond = ggplot(data=df_trial_info,
aes(x=round, y=num_iconic_gestures, fill=condition)) +
geom_boxplot(outlier.shape = NA,
alpha = 0.7) +
scale_fill_manual(values = c("#ED6B06", "#00786A", "darkgrey")) +
scale_y_continuous(limits = c(0, 15), breaks = seq(0, 15, 5)) +
labs(x = "Round",
y = "Mean N of iconic gestures per trial") +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "top",
strip.text = element_text(size = 15, face = 'bold'),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines"))
bp_mean_iconic_by_round_cond
The figure shows that the decline in iconic gestures across rounds does not follow a linear pattern. This suggests that log-transformed round may improve the model fit.
5. [Chosen] Rate of iconic gestures per 100 words
5.1 DataViz: rate per 100 words
bp: mean by condition
bp_mean_iconic_rate_by_cond = ggplot(data=trial_info_pair,
aes(x=condition, y=mean_iconic_per_100words, fill=condition)) +
geom_jitter(aes(color = pair),
size = 1, alpha = 1,
width = 0.07, height = 0) +
geom_boxplot(width = .2,
outlier.shape = NA, alpha = 0.7) +
geom_point(data = trial_info_cond,
aes(y = mean_iconic_per_100words),
shape = 21, size = 3, fill = "white") +
scale_y_continuous(limits = c(0, 8.1), breaks = seq(0, 8, 2)) +
scale_fill_manual(values = c("#ED6B06", "#00786A", "darkgrey")) +
labs(x = "Visibility",
y = "Mean rate of iconic gestures") +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "none",
strip.text = element_text(size = 15, face = 'bold'),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines"))
ggplotly(bp_mean_iconic_rate_by_cond)bp: mean by condition x round
pd = position_dodge(width = .75)
ggplot(data=df_trial_info,
aes(x=round, y=n_iconic_per_100words, fill=condition)) +
geom_boxplot(outlier.shape = NA,
alpha = 0.7) +
geom_point(data = trial_info_cond_round,
aes(y = mean_iconic_per_100words,
group = condition),
position = pd,
shape = 21, size = 2, fill = "white") +
scale_fill_manual(values = c("#ED6B06", "#00786A", "darkgrey")) +
scale_y_continuous(limits = c(0, 22), breaks = seq(0, 20, 5)) +
labs(x = "Round",
y = "Mean rate of iconic gestures") +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "top",
strip.text = element_text(size = 15, face = 'bold'),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines"))
5.2 [Chosen] ZI negative binomial regression models
5.2.1 Prior specification
As the unit of the dependent variable is different from the previous model, we will set different priors for the rate of iconic gestures per 100 words. We will use weakly informative priors for the regression each parameter. To specify priors for the intercept, we will refer to the number of iconic gestures and words reported in Akamine et al. (2024). In the paper, the authors analyzed data from 19 dyads involving in the same task as the current study but in co-present interaction. The total number of iconic gestures reported was 4413, which was collected from 19 dyads, each performing 96 trials. The total number of words produced was 71695 (28152 content) words. Therefore, the expected rate of iconic gestures per 100 words per speaker is \(4413 / (71695 / 100) / 2 = 3.08\) (log(3.08) = 1.12).
For the fixed effects, we will set unbiased priors with a mean of 0 and a standard deviation of 0.5.
For the standard deviation of the random effects, we set the prior to be normal with mean 0 and standard deviation 0.5. For the correlation between the random effects, we set the prior to be LKJ(2).
For the models including the round as fixed effects (whether backward-difference coded or centered + log-transformed), the intercept will represent the mean expected number of iconic gestures (ground mean). As the meaning of the intercept doesn’t change when adding the round variable, we use the same prior for the intercept.
Note that we do not expect the rate of iconic gestures to change across rounds (i.e., we expect the number of words and gestures to decrease at an approximately same pace across the rounds). Also, it is common to set the mean of slopes to be 0 to avoid favoring any directions. Therefore we will set the mean of the prior for slope to 0.
priors_rslope_rate = c(
prior(normal(1.12, 0.5), class = Intercept),
prior(normal(0, 0.5), class = b),
prior(normal(0, 0.5), class = sd),
prior(lkj(2), class = cor))
priors_rslope_rate_zinb = c(
prior(normal(1.12, 0.5), class = Intercept),
prior(normal(0, 0.5), class = b),
prior(normal(0, 0.5), class = sd),
prior(lkj(2), class = cor),
prior(normal(0, 0.5), class = zi), # on the logit scale
prior(normal(0, 50), class = shape))5.2.2 Model 1: condition x round (pair)
In the previous section, we analyzed the number of iconic gestures produced per trial. However, it is common to analyze the rate of iconic gestures per 100 words to account for the differences in the length of the trials and speech rate. Here, we will include the log of speech rate (number of words / 100) as an exposure variable and analyze the rate of iconic gestures per 100 words by condition. Note that the syntax for the exposure variable is different from the Poisson regression model; for negative binomial regression, the exposure variable is included with the rate() function.
Model effect: round
nb_iconic_rate_cond_round = brm(num_iconic_gestures | rate(n_words_100) ~
1 + condition * round +
(1+round|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope_rate,
data = df_trial_info,
sample_prior = T,
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/speakerB/nb_iconic_rate_cond_round")
nb_iconic_rate_cond_round_c = brm(num_iconic_gestures | rate(n_words_100) ~
1 + condition * round_c +
(1+round_c|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope_rate,
data = df_trial_info,
sample_prior = T,
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/speakerB/nb_iconic_rate_cond_round_c")
nb_iconic_rate_cond_round_log = brm(num_iconic_gestures | rate(n_words_100) ~
1 + condition * log_round_c +
(1+log_round_c|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope_rate,
data = df_trial_info,
sample_prior = T,
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/speakerB/nb_iconic_rate_cond_round_log")
### loo compare
if (!file.exists("models/speakerB/loo_nb_iconic_rate_cond_round.rds")){
nb_cond_round_loo = loo(nb_iconic_rate_cond_round)
saveRDS(nb_cond_round_loo, file = "models/speakerB/loo_nb_iconic_rate_cond_round.rds")
}
if (!file.exists("models/speakerB/loo_nb_iconic_rate_cond_round_c.rds")){
nb_cond_round_c_loo = loo(nb_iconic_rate_cond_round_c)
saveRDS(nb_cond_round_c_loo, file = "models/speakerB/loo_nb_iconic_rate_cond_round_c.rds")
}
if (!file.exists("models/speakerB/loo_nb_iconic_rate_cond_round_log.rds")){
nb_cond_round_log_loo = loo(nb_iconic_rate_cond_round_log)
saveRDS(nb_cond_round_log_loo, file = "models/speakerB/loo_nb_iconic_rate_cond_round_log.rds")
}
nb_cond_round_loo = readRDS("models/speakerB/loo_nb_iconic_rate_cond_round.rds")
nb_cond_round_c_loo = readRDS("models/speakerB/loo_nb_iconic_rate_cond_round_c.rds")
nb_cond_round_log_loo = readRDS("models/speakerB/loo_nb_iconic_rate_cond_round_log.rds")
loo_compare(nb_cond_round_loo, nb_cond_round_c_loo, nb_cond_round_log_loo)## elpd_diff se_diff
## nb_iconic_rate_cond_round_c 0.0 0.0
## nb_iconic_rate_cond_round_log -7.8 3.0
## nb_iconic_rate_cond_round -11.5 4.0The leave-one-out (LOO) Effect shows that centered round provides a larger predictive power (elpd_diff) and smaller standard error (se_diff) compared to the backward-difference coded round or centered log-transformed round. Therefore, we will use the centered round.
Model effect: ZI or not
zinb_iconic_rate_cond_round_c = brm(num_iconic_gestures ~ 1 + condition * round_c +
offset(log(n_words_100)) +
(1+round_c|pair) + (1|target),
family = zero_inflated_negbinomial(),
prior = priors_rslope_rate_zinb,
data = df_trial_info,
sample_prior = T,
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/speakerB/zinb_iconic_rate_cond_round_c")
### loo compare
if (!file.exists("models/speakerB/loo_zinb_iconic_rate_cond_round_c.rds")){
zinb_cond_round_c_loo = loo(zinb_iconic_rate_cond_round_c)
saveRDS(zinb_cond_round_c_loo, file = "models/speakerB/loo_zinb_iconic_rate_cond_round_c.rds")
}
nb_cond_round_c_loo = readRDS("models/speakerB/loo_nb_iconic_rate_cond_round_c.rds")
zinb_cond_round_c_loo = readRDS("models/speakerB/loo_zinb_iconic_rate_cond_round_c.rds")
loo_compare(nb_cond_round_c_loo, zinb_cond_round_c_loo)## elpd_diff se_diff
## nb_iconic_rate_cond_round_c 0.0 0.0
## zinb_iconic_rate_cond_round_c -13.1 5.6The leave-one-out (LOO) Effect shows that zero-inflation model has a higher predictive power. As such, we will use zero-inflated negative binomial regression model.
Prior predictive check
zinb_iconic_rate_cond_round_c_prior = brm(num_iconic_gestures ~ 1 + condition * round_c +
offset(log(n_words_100)) +
(1+round_c|pair) + (1|target),
family = zero_inflated_negbinomial(),
prior = priors_rslope_rate_zinb,
sample_prior = "only",
data = df_trial_info,
file = "models/speakerB/zinb_iconic_rate_cond_round_c_prior")
pp_check(zinb_iconic_rate_cond_round_c_prior, ndraws = 100, type = "bars") +
coord_cartesian(xlim = c(0, 20),
ylim = c(0, 4000))
The prior predictive check shows that the model generates data that are somewhat similar to the observed data.
Fit the model
zinb_iconic_rate_cond_round_c = brm(num_iconic_gestures ~ 1 + condition * round_c +
offset(log(n_words_100)) +
(1+round_c|pair) + (1|target),
family = zero_inflated_negbinomial(),
prior = priors_rslope_rate_zinb,
data = df_trial_info,
sample_prior = T,
save_pars = save_pars(all = TRUE),
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/speakerB/zinb_iconic_rate_cond_round_c")
model = zinb_iconic_rate_cond_round_c
summary(model)## Family: zero_inflated_negbinomial
## Links: mu = log; shape = identity; zi = identity
## Formula: num_iconic_gestures ~ 1 + condition * round_c + offset(log(n_words_100)) + (1 + round_c | pair) + (1 | target)
## Data: df_trial_info (Number of observations: 4315)
## Draws: 4 chains, each with iter = 20000; warmup = 2000; thin = 1;
## total post-warmup draws = 72000
##
## Multilevel Hyperparameters:
## ~pair (Number of levels: 45)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 1.07 0.12 0.86 1.32 1.00 22737
## sd(round_c) 0.18 0.03 0.12 0.25 1.00 24819
## cor(Intercept,round_c) 0.64 0.14 0.33 0.85 1.00 45391
## Tail_ESS
## sd(Intercept) 38254
## sd(round_c) 43353
## cor(Intercept,round_c) 52376
##
## ~target (Number of levels: 16)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.09 0.03 0.03 0.16 1.00 20734 21617
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 1.06 0.16 0.74 1.37 1.00 12895
## conditionAO_Asym 0.04 0.30 -0.54 0.64 1.00 17833
## conditionAsym_Sym 0.30 0.30 -0.30 0.88 1.00 16364
## round_c -0.15 0.03 -0.22 -0.08 1.00 26726
## conditionAO_Asym:round_c -0.02 0.07 -0.16 0.12 1.00 27271
## conditionAsym_Sym:round_c -0.05 0.07 -0.19 0.09 1.00 26257
## Tail_ESS
## Intercept 23113
## conditionAO_Asym 29556
## conditionAsym_Sym 28836
## round_c 40828
## conditionAO_Asym:round_c 42495
## conditionAsym_Sym:round_c 41516
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 27.46 14.13 12.22 65.98 1.00 97158 54751
## zi 0.06 0.02 0.03 0.10 1.00 97748 46648
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).bayestestR::hdi(model)# bayestestR::hdi(model, ci = 0.89)The coefficients show that the condition does not have a significant effect on the rate of iconic gestures per 100 words. However, there was a significant decrease in the rate of iconic gestures per 100 words across the rounds. This means that the number of iconic gestures decreased more than the number of words did across the rounds. A formal hypothesis testing will be performed later using Bayes factor.
Visualize the posterior distributions
p = plot_posterior(model, interaction = T)
p
Hypothesis testing: Bayes factor
### varying priors for sensitivity analysis
prior_size = c("xs", "s", "l", "xl")
prior_sd = c(0.1, 0.3, 1, 1.5)
bfs_cond_ao_asym = c()
bfs_cond_asym_sym = c()
bfs_round = c()
for (i in 1:length(prior_sd)){
priors = c(
prior(normal(1.82, 0.5), class = Intercept),
set_prior(paste0("normal(0,", prior_sd[i], ")"), class = "b"),
prior(normal(0, 0.5), class = sd),
prior(lkj(2), class = cor),
prior(normal(0, 0.5), class = zi), # on the logit scale
prior(normal(0, 50), class = shape))
fname = paste0("models/speakerB/zinb_iconic_rate_cond_round_c_", prior_size[i])
fit = brm(num_iconic_gestures ~
1 + condition * round_c +
offset(log(n_words_100)) +
(1+round_c|pair) + (1|target),
family = zero_inflated_negbinomial(),
prior = priors,
data = df_trial_info,
sample_prior = T,
save_pars = save_pars(all = TRUE),
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = fname)
# BF for sym - asym
h = hypothesis(fit, "conditionAO_Asym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_ao_asym = c(bfs_cond_ao_asym, bf)
# BF for sym - ao
h = hypothesis(fit, "conditionAsym_Sym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_asym_sym = c(bfs_cond_asym_sym, bf)
# BF for round
h = hypothesis(fit, "round_c = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round = c(bfs_round, bf)
}
### add BF for the main/medium model
prior_size[3:5] = c("m", prior_size[3:4])
prior_sd[3:5] = c(0.5, prior_sd[3:4])
# BF for sym - asym
h = hypothesis(model, "conditionAO_Asym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_ao_asym[3:5] = c(bf, bfs_cond_ao_asym[3:4])
# BF for sym - ao
h = hypothesis(model, "conditionAsym_Sym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_asym_sym[3:5] = c(bf, bfs_cond_asym_sym[3:4])
# BF for round
h = hypothesis(model, "round_c = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round[3:5] = c(bf, bfs_round[3:4])
### make a df for BFs
df_bf = data.frame(size = prior_size,
sd = prior_sd,
ao_asym = bfs_cond_ao_asym,
asym_sym = bfs_cond_asym_sym,
round = bfs_round) %>%
mutate(prior = paste0("N(0, ", sd, ")")) %>%
pivot_longer(cols = c("ao_asym", "asym_sym", "round"),
names_to = "Effect",
values_to = "BF10") %>%
mutate(Effect = factor(Effect,
levels = c("ao_asym", "asym_sym", "round"),
labels = c("AO-AsymAV", "AsymAV-SymAV", "Round")),
Predictor = ifelse(Effect == "round", "Round", "Visibility"),
BF10_log10 = log10(BF10))
df_bf %>% arrange(Effect, sd)#### Plot BFs ####
ggplot(df_bf, aes(x = prior, y = BF10, group = Effect)) +
geom_hline(yintercept = 1, linetype="dashed") +
geom_point(aes(color=Effect)) +
geom_line(aes(color=Effect)) +
facet_wrap(vars(Predictor), scales="free_y") +
theme_bw(base_size = 12)+
theme(legend.position = "top")+
scale_y_log10("Bayes factor (BF10)",
breaks = c(0.001, 0.03, 0.01, 0.1, 0.33, 1, 3, 10, 30, 100),
labels = c(0.001, 0.03, 0.01, 0.1, 0.33, 1, 3, 10, 30, 100)) +
xlab("prior")
Posterior predictive check
model = zinb_iconic_rate_cond_round_c
pp_check_sym = pp_check_each_condition(model, df_trial_info, "SymAV")
pp_check_asym = pp_check_each_condition(model, df_trial_info, "AsymAV")
pp_check_ao = pp_check_each_condition(model, df_trial_info, "AO")
gridExtra::grid.arrange(pp_check_sym, pp_check_asym, pp_check_ao, ncol = 2)
pp_check_r1 = pp_check_each_round(model, df_trial_info, "R1")
pp_check_r2 = pp_check_each_round(model, df_trial_info, "R2")
pp_check_r3 = pp_check_each_round(model, df_trial_info, "R3")
pp_check_r4 = pp_check_each_round(model, df_trial_info, "R4")
pp_check_r5 = pp_check_each_round(model, df_trial_info, "R5")
pp_check_r6 = pp_check_each_round(model, df_trial_info, "R6")
gridExtra::grid.arrange(pp_check_r1, pp_check_r2, pp_check_r3, pp_check_r4, pp_check_r5, pp_check_r6, ncol = 3)
Pair-wise Effect
emmeans(model, pairwise ~ condition)$contrasts## contrast estimate lower.HPD upper.HPD
## SymAV - AsymAV 0.303 -0.279 0.901
## SymAV - AO 0.347 -0.325 1.041
## AsymAV - AO 0.044 -0.544 0.636
##
## Point estimate displayed: median
## Results are given on the log (not the response) scale.
## HPD interval probability: 0.95emmeans(model, pairwise ~ condition, level = 0.89)$contrasts## contrast estimate lower.HPD upper.HPD
## SymAV - AsymAV 0.303 -0.177 0.785
## SymAV - AO 0.347 -0.215 0.889
## AsymAV - AO 0.044 -0.427 0.534
##
## Point estimate displayed: median
## Results are given on the log (not the response) scale.
## HPD interval probability: 0.89==== Gestural alignment ====
6. Causal model of gestural alignment
Experts in the field of statistics and causal inference have advised that researchers should build a causal model and examine which factors should be included and excluded from regression models if their aim is to infer the causal effects (e.g., McElreath, 2020, Pearl, 2010, Cinelli, Forney, & Pearl, 2020). Following this advice, we will build a causal model to infer the direct effect of speaker visibility on gestural alignment.
We assume the following causal model:
### causal model for gestural alignment rate
dag_gest <- dagitty('dag {
visibility -> gest_align
visibility -> n_words
visibility -> n_gestures
n_words -> lex_align
n_words -> n_gestures
round -> n_words
round -> lex_align
round -> gest_align
round -> n_gestures
role -> n_words
role -> lex_align
role -> gest_align
role -> n_gestures
n_gestures -> gest_align
lex_align -> gest_align
}')
plot(dag_gest)
### check the minimam adjustment set
print("Direct effect of visibility on gestural alignment frequency")## [1] "Direct effect of visibility on gestural alignment frequency"adjustmentSets(dag_gest, exposure = "visibility", outcome = "gest_align", effect = "direct")## { lex_align, n_gestures, role, round }
## { n_gestures, n_words, role, round }print("Direct effect of lexical alignment rate on gestural alignment frequency")## [1] "Direct effect of lexical alignment rate on gestural alignment frequency"adjustmentSets(dag_gest, exposure = "lex_align", outcome = "gest_align", effect = "direct")## { n_gestures, role, round, visibility }
## { n_words, role, round }print("Direct effect of round on gestural alignment frequency")## [1] "Direct effect of round on gestural alignment frequency"adjustmentSets(dag_gest, exposure = "round", outcome = "gest_align", effect = "direct")## { lex_align, n_gestures, role, visibility }The minimum adjustment set for estimating the direct causal effect of speaker visibility on gestural alignment frequency is {lex_align, n_words, round}. Note that because dagitty didn’t find any adjustment set for the direct effect of visibility on lexical alignment frequency when we expected bidirectional causation between lexical and gestural alignment, we assumed a unidirectional causation from lexical alignment to gestural alignment only in this model. This will be reversed in the causal model for lexical alignment frequency, such that we assume a unidirectional causation from gestural alignment to lexical alignment.
In addition, we are also interested in the direct effect of lexical alignment frequency on gestural alignment frequency. The minimum adjustment set for is {visibility, n_gestures, round}.
As the minimum adjustment sets for the direct effects of visibility, lexical alignment frequency, and round on gestural alignment frequency are identical (i.e., {visibility, lex_align, n_gestures, round}), we can estimate the direct effect of these variables on gestural alignment frequency with the following model:
\[ gest\_align \: | \: \text{rate}(n\_iconic\_gesture) \sim \\ \text{visibility} * \text{round} + \text{n_lexical_alignment} + \\ (1+\text{round} | \text{participant}) + (1 | \text{item}) \]
7. Number of gestural alignment
7.1 DataViz: number of gestural alignment
bp: mean by condition
bp_mean_gest_alignment_by_cond = ggplot(data=trial_info_pair,
aes(x=condition, y=mean_num_gest_align, fill=condition)) +
geom_jitter(aes(color = pair),
size = 1, alpha = 1,
width = 0.07, height = 0) +
geom_boxplot(width = .4,
outlier.shape = NA, alpha = 0.7) +
geom_point(data = trial_info_cond,
aes(y = mean_num_gest_align),
size = 2, shape = 4) +
scale_fill_manual(values = c("#ED6B06", "#00786A", "darkgrey")) +
labs(x = "Visibility",
y = "Mean gestural alignment count") +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "none",
strip.text = element_text(size = 15, face = 'bold'),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines"))
ggplotly(bp_mean_gest_alignment_by_cond)7.2 Negative binomial regression
7.2.1 Prior specification
We will set priors for the intercept based on the expected number of gestural alignment reported in Akamine et al. (2024). In the paper, the authors analyzed data from 19 dyads involving in the same task as the current study but in co-present interaction. The total number of gestural alignment reported was 1086, which was collected from 19 dyads, each performing 96 trials. This leads to the expected number of gestural alignment per speaker to be 1086 / (19 * 96) / 2 = 0.3 (log(0.3) = -1.2). As such, we will set the prior for the intercept to be normal with a mean of -1.2 and a standard deviation of 0.5.
For the fixed effects, we will set unbiased priors with a mean of 0. We set different SDs for each effect because the scale of each predictor is different. For N. lexical alignment and N. iconic gestures, we set a prior of Normal(0, 0.2). This means that if the mean N. lex align increases by 1, we expect the mean N. gestural alignment to change by -0.2 to 0.5, with the most likely size of change to be 0. As for the other predictors, we set a prior of Normal(0. 0.5).
We will also specify a prior for the shape parameter to prevent the model from returning divergent transitions. We will set the prior to be normal with a mean of 0 and a wide standard deviation of 50.
For the standard deviation of the random effects, we set the prior to be normal with mean 0 and standard deviation 0.5. For the correlation between the random effects, we set the prior to be LKJ(2).
### pair
priors_rslope_gest_align_zinb = c(
prior(normal(-1.2, 0.5), class = Intercept),
prior(normal(0, 0.5), class = b),
prior(normal(0, 0.2), class = b, coef = "n_iconic_c"),
prior(normal(0, 0.2), class = b, coef = "lex_align_c"),
prior(normal(0, 0.2), class = b, coef = "lex_align_c:conditionAO_Asym"),
prior(normal(0, 0.2), class = b, coef = "lex_align_c:conditionAsym_Sym"),
prior(normal(0, 0.5), class = sd),
prior(lkj(2), class = cor),
prior(normal(0, 0.5), class = zi), # on the logit scale
prior(normal(0, 50), class = shape))7.2.2 Model 2: [ppB] effect of lex align on gest align
Prior predictive check
zinb_gest_align_prior = brm(num_gestural_align ~
1 + lex_align_c * condition + round + n_iconic_c + role +
(1+round|pair) + (1|target),
family = zero_inflated_negbinomial(),
prior = priors_rslope_gest_align_zinb,
data = df_trial_info,
sample_prior = "only",
control = list(adapt_delta = 0.9,
max_treedepth = 20),
file = "models/speakerB/zinb_gest_align_prior")
pp_check(zinb_gest_align_prior, ndraws = 100, type = "bars") +
coord_cartesian(xlim = c(0, 20),
ylim = c(0, 4000))
The prior predictive check shows that the model expects fewer amount of 2, 3, and 4. This suggests that the zero-inflation prior may be too large or the mean for the intercept prior is too low. We will check the prior-posterior update plot and posterior predictive check to see if the model generates data that are similar to the observed data. If not, we will consider modifying the priors.
Fit the model
zinb_align_cond_round = brm(num_gestural_align ~
1 + lex_align_c * condition + round + n_iconic_c + role +
(1+round|pair) + (1|target),
family = zero_inflated_negbinomial(),
prior = priors_rslope_gest_align_zinb,
data = df_trial_info,
sample_prior = T,
save_pars = save_pars(all = TRUE),
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/speakerB/zinb_align_cond_round")
model = zinb_align_cond_round
summary(model)## Family: zero_inflated_negbinomial
## Links: mu = log; shape = identity; zi = identity
## Formula: num_gestural_align ~ 1 + lex_align_c * condition + round + n_iconic_c + role + (1 + round | pair) + (1 | target)
## Data: df_trial_info (Number of observations: 4315)
## Draws: 4 chains, each with iter = 20000; warmup = 2000; thin = 1;
## total post-warmup draws = 72000
##
## Multilevel Hyperparameters:
## ~pair (Number of levels: 45)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 1.13 0.14 0.88 1.44 1.00 17131
## sd(roundR12) 0.43 0.22 0.03 0.86 1.00 17098
## sd(roundR23) 0.23 0.14 0.01 0.52 1.00 19239
## sd(roundR34) 0.29 0.15 0.02 0.60 1.00 16752
## sd(roundR45) 0.15 0.11 0.01 0.41 1.00 32604
## sd(roundR56) 0.31 0.20 0.02 0.73 1.00 20841
## cor(Intercept,roundR12) 0.11 0.28 -0.46 0.63 1.00 64233
## cor(Intercept,roundR23) 0.23 0.31 -0.43 0.75 1.00 54697
## cor(roundR12,roundR23) -0.07 0.32 -0.66 0.58 1.00 51892
## cor(Intercept,roundR34) 0.20 0.29 -0.41 0.71 1.00 59985
## cor(roundR12,roundR34) 0.04 0.31 -0.57 0.62 1.00 40953
## cor(roundR23,roundR34) 0.05 0.32 -0.57 0.64 1.00 37496
## cor(Intercept,roundR45) 0.04 0.33 -0.60 0.66 1.00 85227
## cor(roundR12,roundR45) -0.01 0.33 -0.63 0.62 1.00 76672
## cor(roundR23,roundR45) 0.05 0.33 -0.60 0.67 1.00 62965
## cor(roundR34,roundR45) 0.00 0.33 -0.62 0.63 1.00 70318
## cor(Intercept,roundR56) -0.05 0.31 -0.63 0.55 1.00 77290
## cor(roundR12,roundR56) 0.02 0.32 -0.59 0.63 1.00 58151
## cor(roundR23,roundR56) -0.08 0.33 -0.67 0.57 1.00 51469
## cor(roundR34,roundR56) 0.04 0.32 -0.59 0.64 1.00 55702
## cor(roundR45,roundR56) -0.06 0.34 -0.67 0.60 1.00 44166
## Tail_ESS
## sd(Intercept) 34547
## sd(roundR12) 18949
## sd(roundR23) 22418
## sd(roundR34) 22663
## sd(roundR45) 31508
## sd(roundR56) 26062
## cor(Intercept,roundR12) 52657
## cor(Intercept,roundR23) 51778
## cor(roundR12,roundR23) 53798
## cor(Intercept,roundR34) 47622
## cor(roundR12,roundR34) 50156
## cor(roundR23,roundR34) 52277
## cor(Intercept,roundR45) 54166
## cor(roundR12,roundR45) 55502
## cor(roundR23,roundR45) 57530
## cor(roundR34,roundR45) 59049
## cor(Intercept,roundR56) 52130
## cor(roundR12,roundR56) 55328
## cor(roundR23,roundR56) 53880
## cor(roundR34,roundR56) 58743
## cor(roundR45,roundR56) 55101
##
## ~target (Number of levels: 16)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.20 0.07 0.09 0.35 1.00 22344 23062
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept -2.49 0.18 -2.85 -2.14 1.00
## lex_align_c 0.29 0.04 0.22 0.36 1.00
## conditionAO_Asym 0.38 0.32 -0.26 1.02 1.00
## conditionAsym_Sym 0.16 0.31 -0.45 0.78 1.00
## roundR12 0.97 0.17 0.64 1.31 1.00
## roundR23 -0.04 0.12 -0.29 0.20 1.00
## roundR34 -0.32 0.14 -0.61 -0.05 1.00
## roundR45 -0.30 0.14 -0.57 -0.03 1.00
## roundR56 -0.21 0.17 -0.56 0.13 1.00
## n_iconic_c 0.19 0.02 0.15 0.23 1.00
## role1 -0.99 0.11 -1.20 -0.77 1.00
## lex_align_c:conditionAO_Asym 0.04 0.08 -0.11 0.20 1.00
## lex_align_c:conditionAsym_Sym -0.02 0.07 -0.15 0.11 1.00
## Bulk_ESS Tail_ESS
## Intercept 11238 22288
## lex_align_c 69125 56747
## conditionAO_Asym 17516 31821
## conditionAsym_Sym 15585 27493
## roundR12 59576 53288
## roundR23 50720 49383
## roundR34 46996 52154
## roundR45 65405 55986
## roundR56 65585 52540
## n_iconic_c 35143 37328
## role1 64897 55097
## lex_align_c:conditionAO_Asym 74762 56984
## lex_align_c:conditionAsym_Sym 75600 55472
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 4.57 3.84 2.27 10.66 1.00 31530 24234
## zi 0.04 0.02 0.01 0.09 1.00 65488 42209
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).bayestestR::hdi(model)The coefficients show that the number of lexical alignment is a reliable predictor of the number of gestural alignment.
Hypothesis testing: Bayes factor
### varying priors for sensitivity analysis
prior_size = c("xs", "s", "l", "xl")
prior_sd = c(0.05, 0.1, 0.3, 0.5)
bfs_lex_align = c()
bfs_ao_asym_lex = c()
bfs_asym_sym_lex = c()
for (i in 1:length(prior_sd)){
priors = c(
prior(normal(-1.2, 0.5), class = Intercept),
prior(normal(0, 0.5), class = b),
prior(normal(0, 0.2), class = b, coef = "n_iconic_c"),
prior(normal(0, 0.2), class = b, coef = "lex_align_c"),
prior(normal(0, 0.2), class = b, coef = "lex_align_c:conditionAO_Asym"),
prior(normal(0, 0.2), class = b, coef = "lex_align_c:conditionAsym_Sym"),
prior(normal(0, 0.5), class = sd),
prior(lkj(2), class = cor),
prior(normal(0, 0.5), class = zi),
prior(normal(0, 50), class = shape))
fname = paste0("models/speakerB/zinb_align_cond_round_", prior_size[i])
fit = brm(num_gestural_align ~
1 + lex_align_c * condition + round + n_iconic_c + role +
(1+round|pair) + (1|target),
family = zero_inflated_negbinomial(),
prior = priors,
data = df_trial_info,
sample_prior = T,
save_pars = save_pars(all = TRUE),
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = fname)
### BF for N. lex alignment
h = hypothesis(fit, "lex_align_c = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_lex_align = c(bfs_lex_align, bf)
### BF for interaction
h = hypothesis(fit, "lex_align_c:conditionAO_Asym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_ao_asym_lex = c(bfs_ao_asym_lex, bf)
h = hypothesis(fit, "lex_align_c:conditionAsym_Sym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_asym_sym_lex = c(bfs_asym_sym_lex, bf)
}
### add BF for the main/medium model
prior_size[3:5] = c("m", prior_size[3:4])
prior_sd[3:5] = c(0.2, prior_sd[3:4])
### BF for N. lex alignment
h = hypothesis(model, "lex_align_c = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_lex_align[3:5] = c(bf, bfs_lex_align[3:4])
### BF for interaction
h = hypothesis(model, "lex_align_c:conditionAO_Asym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_ao_asym_lex[3:5] = c(bf, bfs_ao_asym_lex[3:4])
h = hypothesis(model, "lex_align_c:conditionAsym_Sym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_asym_sym_lex[3:5] = c(bf, bfs_asym_sym_lex[3:4])
### make a df for BFs
df_bf_lex = data.frame(size = prior_size,
sd = prior_sd,
lex_align = bfs_lex_align,
ao_asym_lex = bfs_ao_asym_lex,
asym_sym_lex = bfs_asym_sym_lex) %>%
mutate(prior = paste0("N(0, ", sd, ")")) %>%
pivot_longer(cols = c("lex_align", "ao_asym_lex", "asym_sym_lex"),
names_to = "Effect",
values_to = "BF10") %>%
mutate(Effect = ifelse(Effect == "lex_align", "N. lex align",
ifelse(Effect == "ao_asym_lex", "AO--AsymAV:Lex align",
"AsymAV--SymAV: Lex align")),
Predictor = "N. lex align")#### Plot BFs ####
ggplot(filter(df_bf_lex, Effect != "N. lex align"),
aes(x = factor(sd), y = BF10, group = Effect)) +
geom_hline(yintercept = 1, linetype="dashed") +
geom_point(aes(color=Effect)) +
geom_line(aes(color=Effect)) +
facet_wrap(vars(Predictor)) +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "top",
strip.text = element_text(size = 15, face = 'bold'),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines")) +
scale_y_log10("Bayes factor (BF10)",
breaks = c(0.001, 0.03, 0.01, 0.1, 0.33, 1, 3, 10, 30, 100),
labels = c(0.001, 0.03, 0.01, 0.1, 0.33, 1, 3, 10, 30, 100)) +
xlab("SD for the prior")
8. Gestural alignment rate per iconic gestures (n_gest_alignment / n_iconic)
We can analyze the proportion of gestural alignment in two ways: (1) modeling the rate of gestural alignment per iconic gesture using a negative binomial regression model and (2) modeling the probability of gestural alignment using a zero-one-inflated beta regression model.
Key differences in the two models are that the negative binomial regression model assumes that the rate of gestural alignment is a count variable, while the zero-one-inflated beta regression model assumes that the proportion of gestural alignment is a continuous variable bounded between 0 and 1. Also, while negative binomial regression models the rate of events, the zero-one-inflated beta regression models the probability of events. In the case of the proportion of gestural alignment, two models should yield similar results, but it is important to note that the two models are conceptually different.
As it is a common practice to analyze the frequency measures (e.g., number of gestures) as rates, we will use negative binomial regressions.
8.1 DataViz: proportion of gestural alignment
bp: mean by condition
bp_mean_gest_alignment_prop_by_cond = ggplot(data=trial_info_pair,
aes(x=condition, y=mean_gest_align_prop, fill=condition)) +
geom_jitter(aes(color = pair),
size = 1, alpha = 1,
width = 0.07, height = 0) +
geom_boxplot(width = .3,
outlier.shape = NA, alpha = 0.7) +
geom_point(data = trial_info_cond,
aes(y = mean_gest_align_prop),
shape = 21, size = 3, fill = "white") +
scale_fill_manual(values = c("#ED6B06", "#00786A", "darkgrey")) +
labs(x = "Visibility",
y = "Mean gest align rate") +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "none",
strip.text = element_text(size = 15, face = 'bold'),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines"))
ggplotly(bp_mean_gest_alignment_prop_by_cond)bp: mean by condition and round
pd = position_dodge(width = .75)
bp_mean_gest_alignment_prop_by_round_cond =
ggplot(data=trial_info_pair_round,
aes(x=round_n, y=mean_gest_align_prop, fill = condition)) +
geom_jitter(aes(color = pair),
size = 0.5, alpha = 0.7,
width = 0.07, height = 0) +
geom_boxplot(width = .5,
outlier.shape = NA, alpha = 0.7) +
geom_point(data = trial_info_cond_round,
aes(y = mean_gest_align_prop),
position = pd,
shape = 21, size = 2, fill = "white") +
scale_fill_manual(values = c("#ED6B06", "#00786A", "darkgrey")) +
scale_y_continuous(limits = c(0, 1)) +
labs(x = "Round",
y = "Mean gest align rate") +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "none",
strip.text = element_text(size = 15, face = 'bold'),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines")) +
facet_grid(cols = vars(condition))
ggplotly(bp_mean_gest_alignment_prop_by_round_cond)8.2 Prepare data
To model the proportion of gestural alignment, we need to remove trials where the number of iconic gestures is 0. This is because dividing any numbers by 0 results in an undefined value (NA) and prevents model from running.
df_gest_align_posreg_prop = df_trial_info %>%
filter(num_iconic_gestures > 0)
print(paste0("Number of rows before removing trials with no iconic gestures: ", nrow(df_trial_info)))## [1] "Number of rows before removing trials with no iconic gestures: 4315"print(paste0("Number of rows before after trials with no iconic gestures: ", nrow(df_gest_align_posreg_prop)))## [1] "Number of rows before after trials with no iconic gestures: 1264"print(paste0("Number of removed trials: ", nrow(df_trial_info) - nrow(df_gest_align_posreg_prop)))## [1] "Number of removed trials: 3051"8.3 Prior specification
We will set priors based on Akamine et al. (2024). As mentioned in the previous analysis, they detected 4413 iconic gestures and 1086 instances of gestural alignment. Dividing the number of gestural alignment by the number of iconic gestures gives 0.25 (1086/4413) (-1.39 on the log scale). This means that we expect 1 gestural alignment per 4 iconic gestures.
For the slope parameters, we set the mean of 0 with a SD of 0.5. This means that for example, we expect the mean difference between the AO and AsymAV conditions to range from -0.16 to 0.43.
### pair
priors_rint_gest_align_prop = c(
prior(normal(-1.39, 0.5), class = Intercept),
prior(normal(0, 0.5), class = b),
prior(normal(0, 50), class = shape),
prior(normal(0, 0.5), class = sd))
priors_rslope_gest_align_prop = c(
prior(normal(-1.39, 0.5), class = Intercept),
prior(normal(0, 0.5), class = b),
prior(normal(0, 50), class = shape),
prior(normal(0, 0.5), class = sd),
prior(lkj(2), class = cor))
priors_rslope_gest_align_prop_zinb = c(
prior(normal(-1.39, 0.5), class = Intercept),
prior(normal(0, 0.5), class = b),
prior(normal(0, 0.5), class = sd),
prior(lkj(2), class = cor),
prior(normal(0, 0.5), class = zi), # on the logit scale
prior(normal(0, 50), class = shape))8.4 Negative binomial regression models
8.4.1 Model 3: [ppB] condition * round + lex_align_c
Model effect: round
nb_align_rate_cond_round = brm(num_gestural_align | rate(num_iconic_gestures) ~
1 + condition * round + lex_align_c + role +
(1+round|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope_gest_align_prop,
data = df_gest_align_posreg_prop,
sample_prior = T,
save_pars = save_pars(all = TRUE),
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/speakerB/nb_align_rate_cond_round")
nb_align_rate_cond_round_c = brm(num_gestural_align | rate(num_iconic_gestures) ~
1 + condition * round_c + lex_align_c + role +
(1+round_c|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope_gest_align_prop,
data = df_gest_align_posreg_prop,
sample_prior = T,
save_pars = save_pars(all = TRUE),
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/speakerB/nb_align_rate_cond_round_c")
nb_align_rate_cond_round_log = brm(num_gestural_align | rate(num_iconic_gestures) ~
1 + condition * log_round_c + lex_align_c + role +
(1+log_round_c|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope_gest_align_prop,
data = df_gest_align_posreg_prop,
sample_prior = T,
save_pars = save_pars(all = TRUE),
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/speakerB/nb_align_rate_cond_round_log")
### loo compare
if (!file.exists("models/speakerB/loo_nb_align_rate_cond_round.rds")){
nb_cond_round_loo = loo(nb_align_rate_cond_round)
saveRDS(nb_cond_round_loo, file = "models/speakerB/loo_nb_align_rate_cond_round.rds")
}
if (!file.exists("models/speakerB/loo_nb_align_rate_cond_round_c.rds")){
nb_cond_round_c_loo = loo(nb_align_rate_cond_round_c)
saveRDS(nb_cond_round_c_loo, file = "models/speakerB/loo_nb_align_rate_cond_round_c.rds")
}
if (!file.exists("models/speakerB/loo_nb_align_rate_cond_round_log.rds")){
nb_cond_round_log_loo = loo(nb_align_rate_cond_round_log)
saveRDS(nb_cond_round_log_loo, file = "models/speakerB/loo_nb_align_rate_cond_round_log.rds")
}
nb_cond_round_loo = readRDS("models/speakerB/loo_nb_align_rate_cond_round.rds")
nb_cond_round_c_loo = readRDS("models/speakerB/loo_nb_align_rate_cond_round_c.rds")
nb_cond_round_log_loo = readRDS("models/speakerB/loo_nb_align_rate_cond_round_log.rds")
loo_compare(nb_cond_round_loo, nb_cond_round_c_loo, nb_cond_round_log_loo)## elpd_diff se_diff
## nb_align_rate_cond_round 0.0 0.0
## nb_align_rate_cond_round_log -72.5 11.9
## nb_align_rate_cond_round_c -133.9 17.1The leave-one-out (LOO) Effect shows that the backward-difference coded round provides a far larger predictive power (elpd_diff) and a far smaller standard error (se_diff) compared to the centered round or centered log-transformed round. Thus, we will use the bd coded round as a predictor for further analyses.
Model effect: ZI or not
zinb_align_rate_cond_round = brm(num_gestural_align ~
1 + condition * round + lex_align_c + role +
offset(log(num_iconic_gestures)) +
(1+round|pair) + (1|target),
family = zero_inflated_negbinomial(),
prior = priors_rslope_gest_align_prop_zinb,
data = df_gest_align_posreg_prop,
sample_prior = T,
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/speakerB/zinb_align_rate_cond_round")
### loo compare
if (!file.exists("models/speakerB/loo_zinb_align_rate_cond_round.rds")){
zinb_cond_round_c_loo = loo(zinb_align_rate_cond_round)
saveRDS(zinb_cond_round_c_loo, file = "models/speakerB/loo_zinb_align_rate_cond_round.rds")
}
nb_cond_round_c_loo = readRDS("models/speakerB/loo_nb_align_rate_cond_round.rds")
zinb_cond_round_c_loo = readRDS("models/speakerB/loo_zinb_align_rate_cond_round.rds")
loo_compare(nb_cond_round_c_loo, zinb_cond_round_c_loo)## elpd_diff se_diff
## nb_align_rate_cond_round 0.0 0.0
## zinb_align_rate_cond_round -1.4 0.4The loo comparsion shows that non-inflation model has a higher predictive power than the zero-inflation model. As such, we will use NB regression models for further analyses.
Prior predictive check
nb_gest_align_prop_prior = brm(num_gestural_align | rate(num_iconic_gestures) ~
1 + condition + round + lex_align_c + role +
(1+round|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope_gest_align_prop,
data = df_gest_align_posreg_prop,
sample_prior = "only",
control = list(adapt_delta = 0.9,
max_treedepth = 20),
file = "models/speakerB/nb_gest_align_prop_prior")
pp_check(nb_gest_align_prop_prior, ndraws = 100, type = "bars") +
coord_cartesian(xlim = c(0, 20),
ylim = c(0, 3000))
The prior predictive check shows that the model expects fewer amount of 2, 3, and 4. This suggests that the zero-inflation prior may be too large or the mean for the intercept prior is too low. We will check the prior-posterior update plot and posterior predictive check to see if the model generates data that are similar to the observed data. If not, we will consider modifying the priors.
Fit the model
nb_align_rate_cond_round = brm(num_gestural_align | rate(num_iconic_gestures) ~
1 + condition * round + lex_align_c + role +
(1+round|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope_gest_align_prop,
data = df_gest_align_posreg_prop,
sample_prior = T,
save_pars = save_pars(all = TRUE),
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/speakerB/nb_align_rate_cond_round")
model = nb_align_rate_cond_round
summary(model)## Family: negbinomial
## Links: mu = log; shape = identity
## Formula: num_gestural_align | rate(num_iconic_gestures) ~ 1 + condition * round + lex_align_c + role + (1 + round | pair) + (1 | target)
## Data: df_gest_align_posreg_prop (Number of observations: 1264)
## Draws: 4 chains, each with iter = 20000; warmup = 2000; thin = 1;
## total post-warmup draws = 72000
##
## Multilevel Hyperparameters:
## ~pair (Number of levels: 45)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.65 0.11 0.46 0.90 1.00 21086
## sd(roundR12) 0.37 0.17 0.04 0.72 1.00 17598
## sd(roundR23) 0.11 0.08 0.00 0.31 1.00 30412
## sd(roundR34) 0.11 0.09 0.00 0.32 1.00 33413
## sd(roundR45) 0.17 0.12 0.01 0.45 1.00 26547
## sd(roundR56) 0.25 0.17 0.01 0.64 1.00 25128
## cor(Intercept,roundR12) 0.11 0.29 -0.47 0.64 1.00 63968
## cor(Intercept,roundR23) -0.01 0.33 -0.63 0.63 1.00 91834
## cor(roundR12,roundR23) -0.06 0.33 -0.66 0.59 1.00 69885
## cor(Intercept,roundR34) -0.07 0.34 -0.69 0.59 1.00 85044
## cor(roundR12,roundR34) -0.06 0.33 -0.67 0.59 1.00 72391
## cor(roundR23,roundR34) -0.04 0.33 -0.66 0.61 1.00 63507
## cor(Intercept,roundR45) 0.02 0.33 -0.61 0.64 1.00 87625
## cor(roundR12,roundR45) 0.03 0.33 -0.60 0.64 1.00 63588
## cor(roundR23,roundR45) 0.00 0.33 -0.63 0.63 1.00 54163
## cor(roundR34,roundR45) -0.03 0.33 -0.65 0.61 1.00 52871
## cor(Intercept,roundR56) -0.07 0.33 -0.67 0.57 1.00 79111
## cor(roundR12,roundR56) 0.00 0.33 -0.62 0.62 1.00 61216
## cor(roundR23,roundR56) -0.02 0.33 -0.65 0.62 1.00 53120
## cor(roundR34,roundR56) -0.00 0.33 -0.63 0.63 1.00 51432
## cor(roundR45,roundR56) -0.02 0.33 -0.65 0.62 1.00 50780
## Tail_ESS
## sd(Intercept) 37738
## sd(roundR12) 14697
## sd(roundR23) 31325
## sd(roundR34) 32444
## sd(roundR45) 32465
## sd(roundR56) 29115
## cor(Intercept,roundR12) 51402
## cor(Intercept,roundR23) 52909
## cor(roundR12,roundR23) 53430
## cor(Intercept,roundR34) 53772
## cor(roundR12,roundR34) 55527
## cor(roundR23,roundR34) 56199
## cor(Intercept,roundR45) 53607
## cor(roundR12,roundR45) 56174
## cor(roundR23,roundR45) 56077
## cor(roundR34,roundR45) 58994
## cor(Intercept,roundR56) 53182
## cor(roundR12,roundR56) 54035
## cor(roundR23,roundR56) 55885
## cor(roundR34,roundR56) 55300
## cor(roundR45,roundR56) 57829
##
## ~target (Number of levels: 16)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.07 0.05 0.00 0.17 1.00 24855 27141
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept -1.11 0.12 -1.36 -0.88 1.00 16807
## conditionAO_Asym 0.51 0.25 0.02 0.98 1.00 18974
## conditionAsym_Sym 0.02 0.23 -0.44 0.49 1.00 15933
## roundR12 1.70 0.14 1.42 1.97 1.00 56734
## roundR23 0.36 0.10 0.18 0.55 1.00 63234
## roundR34 -0.00 0.11 -0.22 0.21 1.00 56816
## roundR45 -0.10 0.13 -0.36 0.16 1.00 59335
## roundR56 -0.15 0.16 -0.47 0.18 1.00 58659
## lex_align_c 0.05 0.03 -0.00 0.11 1.00 80088
## role1 0.76 0.10 0.57 0.95 1.00 95617
## conditionAO_Asym:roundR12 -0.13 0.29 -0.71 0.44 1.00 60011
## conditionAsym_Sym:roundR12 -0.16 0.25 -0.65 0.35 1.00 57999
## conditionAO_Asym:roundR23 -0.02 0.21 -0.43 0.39 1.00 60027
## conditionAsym_Sym:roundR23 0.22 0.19 -0.16 0.60 1.00 61835
## conditionAO_Asym:roundR34 -0.05 0.23 -0.50 0.41 1.00 56057
## conditionAsym_Sym:roundR34 0.06 0.21 -0.35 0.47 1.00 58737
## conditionAO_Asym:roundR45 0.14 0.27 -0.39 0.67 1.00 63346
## conditionAsym_Sym:roundR45 -0.10 0.24 -0.58 0.38 1.00 61993
## conditionAO_Asym:roundR56 0.42 0.32 -0.21 1.05 1.00 72433
## conditionAsym_Sym:roundR56 0.07 0.28 -0.49 0.62 1.00 73622
## Tail_ESS
## Intercept 31473
## conditionAO_Asym 34023
## conditionAsym_Sym 29667
## roundR12 54551
## roundR23 56243
## roundR34 54612
## roundR45 54190
## roundR56 53561
## lex_align_c 55517
## role1 53785
## conditionAO_Asym:roundR12 52917
## conditionAsym_Sym:roundR12 52033
## conditionAO_Asym:roundR23 54364
## conditionAsym_Sym:roundR23 55893
## conditionAO_Asym:roundR34 54240
## conditionAsym_Sym:roundR34 54102
## conditionAO_Asym:roundR45 57025
## conditionAsym_Sym:roundR45 57614
## conditionAO_Asym:roundR56 57051
## conditionAsym_Sym:roundR56 58041
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 67.42 28.92 23.06 133.77 1.00 98480 49788
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).bayestestR::hdi(model)# bayestestR::hdi(model, ci = 0.89)The coefficients show that the proportion of gestural alignment was significantly higher in the SymAV condition than in the AO condition. Also, the rate of gestural alignment significantly increased from R1–R2 and R2–R3 and stabilized afterwards.
Hypothesis testing: Bayes factor
### varying priors for sensitivity analysis
prior_size = c("xs", "s", "l", "xl")
prior_sd = c(0.1, 0.3, 0.7, 1)
bfs_cond_ao_asym = c()
bfs_cond_asym_sym = c()
bfs_round12 = c()
bfs_round23 = c()
bfs_round34 = c()
bfs_round45 = c()
bfs_round56 = c()
bfs_lex_align = c()
bfs_ao_asym_round12 = c()
bfs_ao_asym_round23 = c()
bfs_ao_asym_round34 = c()
bfs_ao_asym_round45 = c()
bfs_ao_asym_round56 = c()
bfs_asym_sym_round12 = c()
bfs_asym_sym_round23 = c()
bfs_asym_sym_round34 = c()
bfs_asym_sym_round45 = c()
bfs_asym_sym_round56 = c()
for (i in 1:length(prior_sd)){
priors = c(
prior(normal(-1.39, 0.5), class = Intercept),
set_prior(paste0("normal(0,", prior_sd[i], ")"), class = "b"),
prior(normal(0, 0.5), class = sd),
prior(lkj(2), class = cor),
prior(normal(0, 50), class = shape))
fname = paste0("models/speakerB/nb_align_rate_cond_round_", prior_size[i])
fit = brm(num_gestural_align | rate(num_iconic_gestures) ~
1 + condition * round + lex_align_c + role +
(1+round|pair) + (1|target),
family = negbinomial(),
prior = priors,
data = df_gest_align_posreg_prop,
sample_prior = T,
save_pars = save_pars(all = TRUE),
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = fname)
### BF for visibility conditions
# ao - asym
h = hypothesis(fit, "conditionAO_Asym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_ao_asym = c(bfs_cond_ao_asym, bf)
# asym - sym
h = hypothesis(fit, "conditionAsym_Sym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_asym_sym = c(bfs_cond_asym_sym, bf)
### BF for rounds
# R1 - R2
h = hypothesis(model, "roundR12 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round12 = c(bfs_round12, bf)
# R2 - R3
h = hypothesis(model, "roundR23 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round23 = c(bfs_round23, bf)
# R3 - R4
h = hypothesis(model, "roundR34 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round34 = c(bfs_round34, bf)
# R4 - R5
h = hypothesis(model, "roundR45 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round45 = c(bfs_round45, bf)
# R5 - R6
h = hypothesis(model, "roundR56 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round56 = c(bfs_round56, bf)
### BF for N. lex alignment
h = hypothesis(fit, "lex_align_c = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_lex_align = c(bfs_lex_align, bf)
### BF for interaction
# ao - asym: R1 - R2
h = hypothesis(model, "conditionAO_Asym:roundR12 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_ao_asym_round12 = c(bfs_ao_asym_round12, bf)
# ao - asym: R2 - R3
h = hypothesis(model, "conditionAO_Asym:roundR23 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_ao_asym_round23 = c(bfs_ao_asym_round23, bf)
# ao - asym: R3 - R4
h = hypothesis(model, "conditionAO_Asym:roundR34 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_ao_asym_round34 = c(bfs_ao_asym_round34, bf)
# ao - asym: R4 - R5
h = hypothesis(model, "conditionAO_Asym:roundR45 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_ao_asym_round45 = c(bfs_ao_asym_round45, bf)
# ao - asym: R5 - R6
h = hypothesis(model, "conditionAO_Asym:roundR56 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_ao_asym_round56 = c(bfs_ao_asym_round56, bf)
# asym - sym: R1 - R2
h = hypothesis(model, "conditionAsym_Sym:roundR12 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_asym_sym_round12 = c(bfs_asym_sym_round12, bf)
# asym - sym: R2 - R3
h = hypothesis(model, "conditionAsym_Sym:roundR23 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_asym_sym_round23 = c(bfs_asym_sym_round23, bf)
# asym - sym: R3 - R4
h = hypothesis(model, "conditionAsym_Sym:roundR34 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_asym_sym_round34 = c(bfs_asym_sym_round34, bf)
# asym - sym: R4 - R5
h = hypothesis(model, "conditionAsym_Sym:roundR45 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_asym_sym_round45 = c(bfs_asym_sym_round45, bf)
# asym - sym: R5 - R6
h = hypothesis(model, "conditionAsym_Sym:roundR56 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_asym_sym_round56 = c(bfs_asym_sym_round56, bf)
}
### add BF for the main/medium model
prior_size[3:5] = c("m", prior_size[3:4])
prior_sd[3:5] = c(0.5, prior_sd[3:4])
### BF for visibility
# ao - asym
h = hypothesis(model, "conditionAO_Asym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_ao_asym[3:5] = c(bf, bfs_cond_ao_asym[3:4])
# asym - sym
h = hypothesis(model, "conditionAsym_Sym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_asym_sym[3:5] = c(bf, bfs_cond_asym_sym[3:4])
### BF for rounds
# R1 - R2
h = hypothesis(model, "roundR12 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round12[3:5] = c(bf, bfs_round12[3:4])
# R2 - R3
h = hypothesis(model, "roundR23 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round23[3:5] = c(bf, bfs_round23[3:4])
# R3 - R4
h = hypothesis(model, "roundR34 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round34[3:5] = c(bf, bfs_round34[3:4])
# R4 - R5
h = hypothesis(model, "roundR45 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round45[3:5] = c(bf, bfs_round45[3:4])
# R5 - R6
h = hypothesis(model, "roundR56 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round56[3:5] = c(bf, bfs_round56[3:4])
### BF for N. lex alignment
h = hypothesis(model, "lex_align_c = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_lex_align[3:5] = c(bf, bfs_lex_align[3:4])
### BF for interaction
# ao - asym: R1 - R2
h = hypothesis(model, "conditionAO_Asym:roundR12 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_ao_asym_round12[3:5] = c(bf, bfs_ao_asym_round12[3:4])
# ao - asym: R2 - R3
h = hypothesis(model, "conditionAO_Asym:roundR23 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_ao_asym_round23[3:5] = c(bf, bfs_ao_asym_round23[3:4])
# ao - asym: R3 - R4
h = hypothesis(model, "conditionAO_Asym:roundR34 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_ao_asym_round34[3:5] = c(bf, bfs_ao_asym_round34[3:4])
# ao - asym: R4 - R5
h = hypothesis(model, "conditionAO_Asym:roundR45 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_ao_asym_round45[3:5] = c(bf, bfs_ao_asym_round45[3:4])
# ao - asym: R5 - R6
h = hypothesis(model, "conditionAO_Asym:roundR56 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_ao_asym_round56[3:5] = c(bf, bfs_ao_asym_round56[3:4])
# asym - sym: R1 - R2
h = hypothesis(model, "conditionAsym_Sym:roundR12 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_asym_sym_round12[3:5] = c(bf, bfs_asym_sym_round12[3:4])
# asym - sym: R2 - R3
h = hypothesis(model, "conditionAsym_Sym:roundR23 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_asym_sym_round23[3:5] = c(bf, bfs_asym_sym_round23[3:4])
# asym - sym: R3 - R4
h = hypothesis(model, "conditionAsym_Sym:roundR34 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_asym_sym_round34[3:5] = c(bf, bfs_asym_sym_round34[3:4])
# asym - sym: R4 - R5
h = hypothesis(model, "conditionAsym_Sym:roundR45 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_asym_sym_round45[3:5] = c(bf, bfs_asym_sym_round45[3:4])
# asym - sym: R5 - R6
h = hypothesis(model, "conditionAsym_Sym:roundR56 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_asym_sym_round56[3:5] = c(bf, bfs_asym_sym_round56[3:4])
### make a df for BFs
df_bf = data.frame(size = prior_size,
sd = prior_sd,
ao_asym = bfs_cond_ao_asym,
asym_sym = bfs_cond_asym_sym,
round12 = bfs_round12,
round23 = bfs_round23,
round34 = bfs_round34,
round45 = bfs_round45,
round56 = bfs_round56,
lex_align = bfs_lex_align,
ao_asym_round12 = bfs_ao_asym_round12,
ao_asym_round23 = bfs_ao_asym_round23,
ao_asym_round34 = bfs_ao_asym_round34,
ao_asym_round45 = bfs_ao_asym_round45,
ao_asym_round56 = bfs_ao_asym_round56,
asym_sym_round12 = bfs_asym_sym_round12,
asym_sym_round23 = bfs_asym_sym_round23,
asym_sym_round34 = bfs_asym_sym_round34,
asym_sym_round45 = bfs_asym_sym_round45,
asym_sym_round56 = bfs_asym_sym_round56) %>%
mutate(prior = paste0("N(0, ", sd, ")")) %>%
pivot_longer(cols = c("ao_asym", "asym_sym",
"round12", "round23", "round34", "round45", "round56",
"lex_align",
"ao_asym_round12", "ao_asym_round23", "ao_asym_round34",
"ao_asym_round45", "ao_asym_round56",
"asym_sym_round12", "asym_sym_round23", "asym_sym_round34",
"asym_sym_round45", "asym_sym_round56"),
names_to = "Effect",
values_to = "BF10") %>%
mutate(Predictor = ifelse(grepl("_round", Effect), "Interaction",
ifelse(grepl("round", Effect), "Round",
ifelse(Effect == "lex_align", "N. lex align",
"Visibility"))))
df_bf$Effect = recode(df_bf$Effect,
ao_asym = "AO--AsymAV",
asym_sym = "AsymAV--SymAV",
round12 = "R1--R2",
round23 = "R2--R3",
round34 = "R3--R4",
round45 = "R4--R5",
round56 = "R5--R6",
lex_align = "N. lex align",
ao_asym_round12 = "AO--AsymAV:R1--R2",
ao_asym_round23 = "AO--AsymAV:R2--R3",
ao_asym_round34 = "AO--AsymAV:R3--R4",
ao_asym_round45 = "AO--AsymAV:R4--R5",
ao_asym_round56 = "AO--AsymAV:R5--R6",
asym_sym_round12 = "AsymAV--SymAV:R1--R2",
asym_sym_round23 = "AsymAV--SymAV:R2--R3",
asym_sym_round34 = "AsymAV--SymAV:R3--R4",
asym_sym_round45 = "AsymAV--SymAV:R4--R5",
asym_sym_round56 = "AsymAV--SymAV:R5--R6")#### Plot BFs ####
# ggplot(filter(df_bf, Effect!="R1--R2"), #exclude R1--R2 because BF is too huge
# aes(x = factor(sd), y = BF10, group = Effect)) +
# geom_hline(yintercept = 1, linetype="dashed") +
# geom_point(aes(color=Effect)) +
# geom_line(aes(color=Effect)) +
# facet_wrap(vars(Predictor)) +
# theme_clean(base_size = 15) +
# theme(axis.text.x = element_text(colour = "black", size = 14),
# axis.text.y = element_text(colour = "black", size = 14),
# axis.title = element_text(size = 15, face = 'bold'),
# axis.title.x = element_text(vjust = -2),
# axis.title.y = element_text(vjust = 2),
# legend.position = "top",
# strip.text = element_text(size = 15, face = 'bold'),
# plot.background = element_blank(),
# plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines")) +
# scale_y_log10("Bayes factor (BF10)",
# breaks = c(0.001, 0.03, 0.01, 0.1, 0.33, 1, 3, 10, 30, 100),
# labels = c(0.001, 0.03, 0.01, 0.1, 0.33, 1, 3, 10, 30, 100)) +
# xlab("SD for the prior")Pair-wise effect
emmeans(model, pairwise ~ condition)$contrasts## contrast estimate lower.HPD upper.HPD
## SymAV - AsymAV 0.020 -0.446 0.473
## SymAV - AO 0.527 0.015 1.049
## AsymAV - AO 0.507 0.028 0.991
##
## Results are averaged over the levels of: round, role
## Point estimate displayed: median
## Results are given on the log (not the response) scale.
## HPD interval probability: 0.95emmeans(model, pairwise ~ condition, level = 0.89)$contrasts## contrast estimate lower.HPD upper.HPD
## SymAV - AsymAV 0.020 -0.352 0.393
## SymAV - AO 0.527 0.107 0.943
## AsymAV - AO 0.507 0.110 0.892
##
## Results are averaged over the levels of: round, role
## Point estimate displayed: median
## Results are given on the log (not the response) scale.
## HPD interval probability: 0.89Prior-posterior update plot
post_sample = as_draws_df(model)
pp_update_plot(post_sample, model_type="nb")
Posterior predictive check
pp_check_overall = pp_check(model, ndraws = 100, type = "bars") +
coord_cartesian(xlim = c(0, 20),
ylim = c(0, 3000))
pp_check_sym = pp_check_each_condition(model, df_gest_align_posreg_prop, "SymAV")
pp_check_asym = pp_check_each_condition(model, df_gest_align_posreg_prop, "AsymAV")
pp_check_ao = pp_check_each_condition(model, df_gest_align_posreg_prop, "AO")
gridExtra::grid.arrange(pp_check_overall, pp_check_sym,
pp_check_asym, pp_check_ao,
ncol = 2)
Although the model prediction is not perfect, this model had a higher predictive power than the negative binomial model. As such, we will use this model.
Back-transform the effects
### visibility effects
# draw samples
samples = as_draws_df(model)
alpha = samples$b_Intercept
beta_ao_asym = samples$b_conditionAO_Asym
beta_asym_sym = samples$b_conditionAsym_Sym
# convert the samples to the original scale
alpha_orig = exp(alpha)
ao_orig = exp(alpha - beta_ao_asym)
sym_orig = exp(alpha + + beta_asym_sym)
beta_ao_asym_orig = ao_orig - alpha_orig
beta_asym_sym_orig = sym_orig - alpha_orig
# summarize the mean and CI of the samples as dataframe
tibble(estimates = c("intercept_asym", "ao", "sym", "beta_ao_asym", "beta_asym_sym"),
mean = c(mean(alpha_orig), mean(ao_orig), mean(sym_orig),
mean(beta_ao_asym_orig), mean(beta_asym_sym_orig)),
ci_low = c(quantile(alpha_orig, 0.025), quantile(ao_orig, 0.025),
quantile(sym_orig, 0.025), quantile(beta_ao_asym_orig, 0.025),
quantile(beta_asym_sym_orig, 0.025)),
ci_high = c(quantile(alpha_orig, 0.975), quantile(ao_orig, 0.975),
quantile(sym_orig, 0.975), quantile(beta_ao_asym_orig, 0.975),
quantile(beta_asym_sym_orig, 0.975)))8.4.2 Model 4: [ppB] condition_sum * round + lex_align_c
Here, we will run another model to compare the rate of gestural alignment between the SymAV and AO conditions.
Prior predictive check
nb_gest_align_prop_sum_prior = brm(num_gestural_align | rate(num_iconic_gestures) ~
1 + condition_sum + round + lex_align_c + role +
(1+round|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope_gest_align_prop,
data = df_gest_align_posreg_prop,
sample_prior = "only",
control = list(adapt_delta = 0.9,
max_treedepth = 20),
file = "models/speakerB/nb_gest_align_prop_sum_prior")
pp_check(nb_gest_align_prop_sum_prior, ndraws = 100, type = "bars") +
coord_cartesian(xlim = c(0, 20),
ylim = c(0, 3000))
The prior predictive check shows that the model expects fewer amount of 2, 3, and 4. This suggests that the zero-inflation prior may be too large or the mean for the intercept prior is too low. We will check the prior-posterior update plot and posterior predictive check to see if the model generates data that are similar to the observed data. If not, we will consider modifying the priors.
Fit the model
nb_align_rate_cond_round_sum = brm(num_gestural_align | rate(num_iconic_gestures) ~
1 + condition_sum * round + lex_align_c + role +
(1+round|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope_gest_align_prop,
data = df_gest_align_posreg_prop,
sample_prior = T,
save_pars = save_pars(all = TRUE),
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/speakerB/nb_align_rate_cond_round_sum")
model = nb_align_rate_cond_round_sum
summary(model)## Family: negbinomial
## Links: mu = log; shape = identity
## Formula: num_gestural_align | rate(num_iconic_gestures) ~ 1 + condition_sum * round + lex_align_c + role + (1 + round | pair) + (1 | target)
## Data: df_gest_align_posreg_prop (Number of observations: 1264)
## Draws: 4 chains, each with iter = 20000; warmup = 2000; thin = 1;
## total post-warmup draws = 72000
##
## Multilevel Hyperparameters:
## ~pair (Number of levels: 45)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.65 0.11 0.46 0.89 1.00 20133
## sd(roundR12) 0.37 0.17 0.04 0.72 1.00 20020
## sd(roundR23) 0.11 0.08 0.00 0.31 1.00 34135
## sd(roundR34) 0.11 0.09 0.00 0.32 1.00 33498
## sd(roundR45) 0.17 0.12 0.01 0.45 1.00 29102
## sd(roundR56) 0.25 0.18 0.01 0.65 1.00 25173
## cor(Intercept,roundR12) 0.11 0.29 -0.47 0.64 1.00 69482
## cor(Intercept,roundR23) -0.01 0.33 -0.64 0.62 1.00 93543
## cor(roundR12,roundR23) -0.06 0.33 -0.66 0.59 1.00 72794
## cor(Intercept,roundR34) -0.07 0.34 -0.69 0.59 1.00 93955
## cor(roundR12,roundR34) -0.05 0.33 -0.66 0.59 1.00 73807
## cor(roundR23,roundR34) -0.03 0.33 -0.66 0.61 1.00 65352
## cor(Intercept,roundR45) 0.02 0.33 -0.61 0.63 1.00 100701
## cor(roundR12,roundR45) 0.03 0.32 -0.60 0.64 1.00 69318
## cor(roundR23,roundR45) 0.00 0.33 -0.62 0.63 1.00 55490
## cor(roundR34,roundR45) -0.03 0.33 -0.65 0.61 1.00 52700
## cor(Intercept,roundR56) -0.08 0.33 -0.68 0.56 1.00 85654
## cor(roundR12,roundR56) -0.01 0.32 -0.62 0.61 1.00 63761
## cor(roundR23,roundR56) -0.02 0.33 -0.64 0.62 1.00 54273
## cor(roundR34,roundR56) -0.01 0.33 -0.64 0.63 1.00 51739
## cor(roundR45,roundR56) -0.02 0.33 -0.64 0.61 1.00 51362
## Tail_ESS
## sd(Intercept) 34892
## sd(roundR12) 19453
## sd(roundR23) 34420
## sd(roundR34) 33202
## sd(roundR45) 32548
## sd(roundR56) 30242
## cor(Intercept,roundR12) 52311
## cor(Intercept,roundR23) 55027
## cor(roundR12,roundR23) 56339
## cor(Intercept,roundR34) 51632
## cor(roundR12,roundR34) 54064
## cor(roundR23,roundR34) 57292
## cor(Intercept,roundR45) 55336
## cor(roundR12,roundR45) 54192
## cor(roundR23,roundR45) 55671
## cor(roundR34,roundR45) 56635
## cor(Intercept,roundR56) 54051
## cor(roundR12,roundR56) 55655
## cor(roundR23,roundR56) 56716
## cor(roundR34,roundR56) 59091
## cor(roundR45,roundR56) 57312
##
## ~target (Number of levels: 16)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.07 0.05 0.00 0.17 1.00 25716 27728
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept -1.11 0.12 -1.36 -0.88 1.00
## condition_sumAO_Sym 0.48 0.24 -0.01 0.95 1.00
## condition_sumAsym_Sym -0.10 0.23 -0.56 0.36 1.00
## roundR12 1.69 0.14 1.41 1.97 1.00
## roundR23 0.37 0.09 0.18 0.55 1.00
## roundR34 -0.00 0.11 -0.22 0.21 1.00
## roundR45 -0.10 0.13 -0.36 0.16 1.00
## roundR56 -0.15 0.17 -0.47 0.18 1.00
## lex_align_c 0.05 0.03 -0.00 0.11 1.00
## role1 0.76 0.10 0.57 0.95 1.00
## condition_sumAO_Sym:roundR12 -0.21 0.29 -0.78 0.37 1.00
## condition_sumAsym_Sym:roundR12 -0.11 0.26 -0.61 0.41 1.00
## condition_sumAO_Sym:roundR23 0.15 0.21 -0.25 0.57 1.00
## condition_sumAsym_Sym:roundR23 0.20 0.19 -0.18 0.57 1.00
## condition_sumAO_Sym:roundR34 0.02 0.23 -0.43 0.46 1.00
## condition_sumAsym_Sym:roundR34 0.08 0.21 -0.34 0.50 1.00
## condition_sumAO_Sym:roundR45 0.08 0.27 -0.44 0.60 1.00
## condition_sumAsym_Sym:roundR45 -0.09 0.24 -0.57 0.39 1.00
## condition_sumAO_Sym:roundR56 0.41 0.33 -0.24 1.05 1.00
## condition_sumAsym_Sym:roundR56 -0.06 0.28 -0.62 0.49 1.00
## Bulk_ESS Tail_ESS
## Intercept 17891 33150
## condition_sumAO_Sym 18128 31916
## condition_sumAsym_Sym 15984 29275
## roundR12 60269 52721
## roundR23 68573 57116
## roundR34 64296 56352
## roundR45 63043 56515
## roundR56 62677 53499
## lex_align_c 81773 56619
## role1 102399 54135
## condition_sumAO_Sym:roundR12 59957 53532
## condition_sumAsym_Sym:roundR12 61536 54591
## condition_sumAO_Sym:roundR23 65983 56053
## condition_sumAsym_Sym:roundR23 65743 55998
## condition_sumAO_Sym:roundR34 61413 55727
## condition_sumAsym_Sym:roundR34 61294 56291
## condition_sumAO_Sym:roundR45 63529 54058
## condition_sumAsym_Sym:roundR45 64047 56071
## condition_sumAO_Sym:roundR56 68951 56265
## condition_sumAsym_Sym:roundR56 74675 54323
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 67.35 28.95 22.99 133.91 1.00 112593 50071
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).bayestestR::hdi(model)# bayestestR::hdi(model, ci = 0.89)The coefficients show that the proportion of gestural alignment was significantly higher in the SymAV condition than in the AO condition. Also, the rate of gestural alignment significantly increased from R1–R2 and R2–R3 and stabilized afterwards.
Visualize posterior distributions
main_model = nb_align_rate_cond_round
model1 = zinb_align_cond_round
plot_posterior(main_model, model1, model)
Hypothesis testing: Bayes factor
### varying priors for sensitivity analysis
prior_size = c("xs", "s", "l", "xl")
prior_sd = c(0.1, 0.3, 0.7, 1)
bfs_cond_ao_sym = c()
bfs_ao_sym_round12 = c()
bfs_ao_sym_round23 = c()
bfs_ao_sym_round34 = c()
bfs_ao_sym_round45 = c()
bfs_ao_sym_round56 = c()
for (i in 1:length(prior_sd)){
priors = c(
prior(normal(-1.39, 0.5), class = Intercept),
set_prior(paste0("normal(0,", prior_sd[i], ")"), class = "b"),
prior(normal(0, 0.5), class = sd),
prior(lkj(2), class = cor),
prior(normal(0, 50), class = shape))
fname = paste0("models/speakerB/nb_align_rate_cond_round_sum_", prior_size[i])
fit = brm(num_gestural_align | rate(num_iconic_gestures) ~
1 + condition_sum * round + lex_align_c + role +
(1+round|pair) + (1|target),
family = negbinomial(),
prior = priors,
data = df_gest_align_posreg_prop,
sample_prior = T,
save_pars = save_pars(all = TRUE),
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = fname)
### BF for visibility conditions
# BF for ao - sym
h = hypothesis(fit, "condition_sumAO_Sym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_ao_sym = c(bfs_cond_ao_sym, bf)
### BF for interaction
# ao - sym: R1 - R2
h = hypothesis(model, "condition_sumAO_Sym:roundR12 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_ao_sym_round12 = c(bfs_ao_sym_round12, bf)
# ao - sym: R2 - R3
h = hypothesis(model, "condition_sumAO_Sym:roundR23 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_ao_sym_round23 = c(bfs_ao_sym_round23, bf)
# ao - sym: R3 - R4
h = hypothesis(model, "condition_sumAO_Sym:roundR34 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_ao_sym_round34 = c(bfs_ao_sym_round34, bf)
# ao - sym: R4 - R5
h = hypothesis(model, "condition_sumAO_Sym:roundR45 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_ao_sym_round45 = c(bfs_ao_sym_round45, bf)
# ao - sym: R5 - R6
h = hypothesis(model, "condition_sumAO_Sym:roundR56 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_ao_sym_round56 = c(bfs_ao_sym_round56, bf)
}
### add BF for the main/medium model
prior_size[3:5] = c("m", prior_size[3:4])
prior_sd[3:5] = c(0.5, prior_sd[3:4])
### BF for visibility conditions
# ao - sym
h = hypothesis(model, "condition_sumAO_Sym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_ao_sym[3:5] = c(bf, bfs_cond_ao_sym[3:4])
### BF for interaction
# ao - sym: R1 - R2
h = hypothesis(model, "condition_sumAO_Sym:roundR12 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_ao_sym_round12[3:5] = c(bf, bfs_ao_sym_round12[3:4])
# ao - sym: R2 - R3
h = hypothesis(model, "condition_sumAO_Sym:roundR23 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_ao_sym_round23[3:5] = c(bf, bfs_ao_sym_round23[3:4])
# ao - sym: R3 - R4
h = hypothesis(model, "condition_sumAO_Sym:roundR34 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_ao_sym_round34[3:5] = c(bf, bfs_ao_sym_round34[3:4])
# ao - sym: R4 - R5
h = hypothesis(model, "condition_sumAO_Sym:roundR45 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_ao_sym_round45[3:5] = c(bf, bfs_ao_sym_round45[3:4])
# ao - sym: R5 - R6
h = hypothesis(model, "condition_sumAO_Sym:roundR56 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_ao_sym_round56[3:5] = c(bf, bfs_ao_sym_round56[3:4])
### make a df for BFs
# df_bf_temp = data.frame(size = prior_size,
# sd = prior_sd,
# Effect = "AO--SymAV",
# BF10 = bfs_cond_ao_sym) %>%
# mutate(prior = paste0("N(0, ", sd, ")"),
# Predictor = "Visibility")
df_bf_temp = data.frame(size = prior_size,
sd = prior_sd,
ao_sym = bfs_cond_ao_sym,
ao_sym_round12 = bfs_ao_sym_round12,
ao_sym_round23 = bfs_ao_sym_round23,
ao_sym_round34 = bfs_ao_sym_round34,
ao_sym_round45 = bfs_ao_sym_round45,
ao_sym_round56 = bfs_ao_sym_round56) %>%
mutate(prior = paste0("N(0, ", sd, ")")) %>%
pivot_longer(cols = c("ao_sym",
"ao_sym_round12", "ao_sym_round23", "ao_sym_round34",
"ao_sym_round45", "ao_sym_round56"),
names_to = "Effect",
values_to = "BF10") %>%
mutate(Predictor = ifelse(grepl("_round", Effect), "Interaction", "Visibility"))
df_bf_temp$Effect = recode(df_bf_temp$Effect,
ao_sym = "AO--SymAV",
ao_sym_round12 = "AO--SymAV:R1--R2",
ao_sym_round23 = "AO--SymAV:R2--R3",
ao_sym_round34 = "AO--SymAV:R3--R4",
ao_sym_round45 = "AO--SymAV:R4--R5",
ao_sym_round56 = "AO--SymAV:R5--R6")
df_bf_new = df_bf %>%
filter(Effect != "N. lex align") %>%
rbind(df_bf_temp) %>%
rbind(df_bf_lex) %>%
filter(!(Predictor == "N. lex align" & Effect != "N. lex align")) %>%
mutate(Effect = factor(Effect,
levels = c("AO--SymAV", "AO--AsymAV", "AsymAV--SymAV",
"R1--R2", "R2--R3", "R3--R4", "R4--R5", "R5--R6",
"N. lex align",
"AO--SymAV:R1--R2", "AO--SymAV:R2--R3", "AO--SymAV:R3--R4",
"AO--SymAV:R4--R5", "AO--SymAV:R5--R6",
"AO--AsymAV:R1--R2", "AO--AsymAV:R2--R3", "AO--AsymAV:R3--R4",
"AO--AsymAV:R4--R5", "AO--AsymAV:R5--R6",
"AsymAV--SymAV:R1--R2", "AsymAV--SymAV:R2--R3", "AsymAV--SymAV:R3--R4",
"AsymAV--SymAV:R4--R5", "AsymAV--SymAV:R5--R6")),
Predictor = factor(Predictor,
levels = c("Visibility", "Round", "N. lex align", "Interaction")))
df_bf_new %>% arrange(Effect, sd)#### Plot BFs ####
ggplot(
filter(df_bf_new,
Effect != "R1--R2", Predictor != "N. lex align", #exclude effects/predictors where BF is too large
Predictor != "Interaction"), #exclude interactions
aes(x = factor(sd), y = BF10, group = Effect)) +
geom_hline(yintercept = 1, linetype="dashed") +
geom_point(aes(color=Effect)) +
geom_line(aes(color=Effect)) +
facet_wrap(vars(Predictor), scales = "free_x") +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
# legend.position = "top",
legend.title=element_text(size=14),
legend.text=element_text(size=13),
strip.text = element_text(size = 15, face = 'bold'),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines")) +
scale_y_log10("Bayes factor (BF10)",
breaks = c(0.001, 0.03, 0.01, 0.1, 0.33, 1, 3, 10, 30, 100),
labels = c(0.001, 0.03, 0.01, 0.1, 0.33, 1, 3, 10, 30, 100)) +
xlab("SD for the prior")
Probability of direction
p_direction(model)Prior-posterior update plot
post_sample = as_draws_df(model)
pp_update_plot(post_sample, model_type="nb")
Posterior predictive check
pp_check_overall = pp_check(model, ndraws = 100, type = "bars") +
coord_cartesian(xlim = c(0, 20),
ylim = c(0, 3000))
pp_check_sym = pp_check_each_condition(model, df_gest_align_posreg_prop, "SymAV")
pp_check_asym = pp_check_each_condition(model, df_gest_align_posreg_prop, "AsymAV")
pp_check_ao = pp_check_each_condition(model, df_gest_align_posreg_prop, "AO")
gridExtra::grid.arrange(pp_check_overall, pp_check_sym,
pp_check_asym, pp_check_ao,
ncol = 2)
Although the model prediction is not perfect, this model had a higher predictive power than the negative binomial model. As such, we will use this model.
9. Correlation btw lexical and gestural alignment
cor = pcor.test(df_trial_info$lex_align_c, df_trial_info$num_gestural_align, df_trial_info$log_round_c)
cor=====Session info=====
sessionInfo()## R version 4.2.2 (2022-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 26200)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] rstan_2.32.6 StanHeaders_2.32.7 svglite_2.2.2 ppcor_1.1
## [5] MASS_7.3-58.1 doParallel_1.0.17 iterators_1.0.14 foreach_1.5.2
## [9] emmeans_1.10.7 tidybayes_3.0.7 bayestestR_0.15.2 brms_2.22.0
## [13] Rcpp_1.0.12 plotly_4.10.4 rsvg_2.6.2 DiagrammeRsvg_0.1
## [17] DiagrammeR_1.0.11 dagitty_0.3-4 ggh4x_0.3.1 ggthemes_5.1.0
## [21] hypr_0.2.8 plotrix_3.8-4 lubridate_1.9.3 forcats_1.0.0
## [25] stringr_1.5.1 dplyr_1.1.4 purrr_1.0.2 readr_2.1.5
## [29] tidyr_1.3.1 tibble_3.2.1 ggplot2_3.5.2 tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] colorspace_2.1-0 estimability_1.5.1 QuickJSR_1.1.3
## [4] rstudioapi_0.17.1 farver_2.1.1 svUnit_1.0.6
## [7] bit64_4.0.5 mvtnorm_1.2-4 bridgesampling_1.1-2
## [10] codetools_0.2-18 cachem_1.0.8 knitr_1.50
## [13] bayesplot_1.11.1 jsonlite_1.8.8 ggdist_3.3.2
## [16] compiler_4.2.2 httr_1.4.7 backports_1.4.1
## [19] Matrix_1.5-1 fastmap_1.1.1 lazyeval_0.2.2
## [22] cli_3.6.2 visNetwork_2.1.2 htmltools_0.5.8.1
## [25] tools_4.2.2 coda_0.19-4.1 gtable_0.3.6
## [28] glue_1.7.0 reshape2_1.4.4 posterior_1.6.1
## [31] V8_6.0.6 jquerylib_0.1.4 vctrs_0.6.5
## [34] nlme_3.1-160 crosstalk_1.2.1 insight_1.1.0
## [37] tensorA_0.36.2.1 xfun_0.53 ps_1.7.6
## [40] timechange_0.3.0 lifecycle_1.0.4 scales_1.3.0
## [43] vroom_1.6.5 ragg_1.3.0 hms_1.1.3
## [46] Brobdingnag_1.2-9 inline_0.3.21 RColorBrewer_1.1-3
## [49] yaml_2.3.8 curl_6.2.1 gridExtra_2.3
## [52] loo_2.8.0 sass_0.4.9 stringi_1.8.3
## [55] checkmate_2.3.1 pkgbuild_1.4.6 boot_1.3-28
## [58] cmdstanr_0.8.1 rlang_1.1.3 pkgconfig_2.0.3
## [61] systemfonts_1.3.1 matrixStats_1.3.0 distributional_0.5.0
## [64] pracma_2.4.4 evaluate_1.0.3 lattice_0.20-45
## [67] rstantools_2.4.0 htmlwidgets_1.6.4 labeling_0.4.3
## [70] processx_3.8.4 bit_4.0.5 tidyselect_1.2.1
## [73] plyr_1.8.9 magrittr_2.0.3 R6_2.6.1
## [76] generics_0.1.3 pillar_1.10.1 withr_3.0.2
## [79] datawizard_1.0.1 abind_1.4-8 crayon_1.5.3
## [82] arrayhelpers_1.1-0 tzdb_0.4.0 rmarkdown_2.29
## [85] grid_4.2.2 data.table_1.15.4 digest_0.6.35
## [88] xtable_1.8-4 textshaping_0.3.7 stats4_4.2.2
## [91] RcppParallel_5.1.7 munsell_0.5.1 viridisLite_0.4.2
## [94] bslib_0.9.0